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ABSTRACT 


This report describes the study and development 
of two numerical techniques for the analysis of 

electromagnetic scattering from a rectangular wire mesh. Both 
techniques follow from one basic formulation and they are 
both solved in the spectral domain. These techniques were 
developed as a result of an investigation towards more 
efficient numerical computation for mesh scattering. These 
techniques are efficient for the following reasons: 

a) They make use of the Fast Fourier Transform. 

b) They avoid any convolution problems by converting 
integrodifferential equations into algebraic 
equations. 

c) They do not require inversions of any matrices. 

The first method, the "SIT" or Spectral Iteration Technique, 
is applied for regions where the spacing between wires is 
not less than two wavelenghs. The second method, the 
"SDGC" or Spectral Domain Conjugate Gradient approach, can 
be used for any spacing between adjacent wires. A study of 
electromagnetic wave properties, such as reflection 
coefficient, induced currents and aperture fields, as 
functions of frequency, angle of incidence, polarization and 
thickness of wires is presented. Examples and comparisons 
of results with other methods are also included to support 
the validity of the new algorithms. 
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INTRODUCTION 


A new technology for large space-based systems requires 
antennas with 100 meters or larger in diameter for radio 
frequency operation, communication, earth observation and 
radio astronomy applications. 

A new type of antenna, the MESH DEPLOYABLE ANTENNA, 
which appears to be more cost-effective and easier to . 
transport into space compared to a solid reflector of 100 
meters in diameter, was the motivation for the study 
reported herein. The mesh used to construct large space 
reflector antennas is usually made of gold-plated molybdenum 
wire about one mill in diameter. The wires run and cross in 
a weave that is periodic in nature, forming a reflecting 
surface that behaves differently depending on the number of 
openings per wavelength and polarization of the incident 
energy. The undesirable effects resulting from such a 
surface include transmission loss, resistive loss, and cross 
polarization loss. 

Here a wire mesh structure is used as a simplified model 
of the knitted (woven) material. A rectangular mesh 
structure is a periodic structure, and scattering from 
periodic structures is a subject that has a long and 
illustrious history dating back to Lamb and Rayleigh in the 
last century [1-5] . 
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Constructing solutions to the problem of mesh scattering 
can be achieved using a variety of methods. One possible 
method which has been widely used is the METHOD OF MOMENTS 
(MOM) [6-8] . This method, when applied to periodic 
surfaces, has the disadvantage of requiring the inversion of 
a very large matrix, a fact that renders the method 
unwieldy. Other methods involve COUPLED INTEGRAL EQUATIONS. 
These methods will not usually yield a solution due to the 
complexity of inverting the integrals for a periodic mesh. 
Another popular technique used for estimating the reflection 
coefficient from a wire mesh is based on AVERAGED BOUNDARY 
CONDITIONS [9-10] . This method offers good results when the 
number of mesh openings per wavelength is large [11] . 
However, even this method fails for certain applications 
when the number of openings per wavelength becomes small. 

This dissertation includes the analysis and formulation 
of two new models for studying scattering from wire meshes 
that are more efficient and simpler to apply than the 
previous methods. The first method is based on the SPECTRAL 
ITERATION APPROACH (SIT) [12-18] which is valid for cases 
where the spacing between adjacent wires is larger than two 
wavelengths. This limitation on the size of spacing between 
wires for the SIT method led to the development of the 
second model which is valid for all spacings. This new 
model is the SPECTRAL DOMAIN CONJUGATE GRADIENT method 
( SDCG ) [19-22] and is a combination of the SIT and the 
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Conjugate Gradient method. Both methods utilize the fast 
Fourier transform and avoid any convolution problems and any 
inversion of matrices. 

These two techniques offer new accurate models which can 
be extended and applied to the more difficult problems of 
knitted mesh surfaces. A number of examples are computed 
and compared with other methods. Also, comments and 
suggestions are made for possible extension of the SDCG 
method to the more complicated problem of the knitted 


structure . 
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2. THE SPECTRAL ITERATION APPROACH 
2.1 FORMULATION 

Any scattering problem could be expressed in the form of 
the integral equation: 

®(x)= |k(x,x’) 'F(x') dx + <6 inc ( x ) (2.1) 

with the constitutive equation ^(x)=K(x) <3>(x) (2.2) 

where K(x,x') is the kernel of the integral transform 

inc 

® (x) is the externally applied field 

® (x) is the field quantity, and 

¥(x) is the source quantity 

The S.I.T. method is a frequency domain (Spectral 
Domain) solution, and consists of casting the general basic 
global equations (i.e. the second order partial differential 
equation or its integral representation, such as equation 

(2.1) ) as a local algebraic equation in the Fourier 
transform space and leaving the local constitutive equation 
as a local algebraic equation in real space. That is, 
taking the Fourier transform of equation (2.1) and keeping 

(2.2) the same, one arrives at: 

®(k) = K ( k ) W(k) + 5 inc (k) (2.3) 

W (x) = K(x) © (x) (2.4) 

Equations (2.3) and (2.4) show how the original set of 
equations are converted into a set of two simultaneous 
algebraic equations in two unknowns (the fields and the 
induced currents) in two different domains connected by the 
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Fourier transform which is given by: 

- rcO A A 

F(k) = I f ( x ) exp(jk.x) dx (2.5) 

■too 

The operation in equation (2.5) from now on will be denoted 
by the transform pair: 

FU)-*— ^f(x) (2.6) 

By virtue of the numerical Fast Fourier transform and the 
local algebraic representation, the number of required 
complex multiply and add operations and the number of 
required storage locations are of the order of Nlog 2 N and N 
respectively (where N represents the number of Floquet modes 
or cells into which the problem is discretized). 

For periodic structures, the Floquet theorem [23] is 
used to account for the periodicity of the wire mesh and the 
coupling between adjacent wires. The specific equations for 
a wire mesh (See Figure 2.1) are formulated as follows: 



Fig. 2.1. Wire mesh geometry 
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The electric field E due to a magnetic current M is 
given by: 

E(x,y)=-l/ e VxF(x,y,z) (2.7) 

where F* is the associated electric vector potential of the 
source and e is the permitivity of the medium in which the 
wire mesh is placed. T and ift are related by the free space 
Green's function G=exp(-jk.r)/4 ti r by 

"f(T)» Jg( r\ T' ) M(?') dr' (2.8) 

From this the, magnetic field intensity, ft, can be derived 
(See Appendix 8.1) from Maxwell's equations and expressed 
as: 

H ( x ,y , z ) = -ja)F(x,y ,z) + Vy.T(x,y,z)/;j<ijeu (2.9) 

where u is the permeability of the medium. Since we have a 

planar structure F is set equal to zero. Now expanding 

z 

equation (2.9) in terms of its Cartesian coordinates x and y 
yields: 


H( x ,y ) =- 


3ajy.e 


2 

(k_ + 


2 2 

a . a 


Cx dy dx 2 


)F x x +(k + 


2 2 

a a 

° dxdSdy"* 


( 2 . 10 ) 


A planar periodic structure such as that shown in Figure 
(2.1) could be considered to be the source distribution for 
the magnetic field of the equation (2.10). Substituting 
equation (2.8) into equation (2.10) and taking Fourier 
transform of equation (2.10) yields the transformed 
scattered tangential fields at z=0 in the following form: 


^<> 



DCf* 
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3 mn = 2 tx n/c - 2 it m/a cotfl -ko sinO sincp are the Floquet 
modes [24] and 

G ( a , 0 ) = - j/2 ( k 2 - a 2 - 0 2 mri ) I is the Fourier 

mn mn J o mn mn 

transform of Green's function. 



Fig. 2.2 FSS Geometry 

Taking the inverse Fourier transform of equation (2.11) 
yields : 
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n ( X ,y ) = 

CO j U ( 


mn 


k 2 - a 2 a 0 

o mn mn mn 


-a 0 k 2 - 0 2 

mn mn o mn 


G( a , 0 ) M ( a__ , 0 ) 

mn mn mn mn 


.exp [ j ( cl.x + 0„„y)] 
mn mn 


(2.12-a) 

Now, by using the equivalence theorem and applying the 
appropriate boundary conditions on H s (x,y) at 2=0 (See 
Appendix 8.2) leads to: 



. 2 2 
k - a. 
o mn 


~ a mn ^ mnl 


G ( a , 0 __ ) E( a , 0 __ ) 
mn mn mn mn 


• exp [ j ( a ra x + 0 __y) ] 
mn mn- 1 


(2.12-b) 


where E represents the transformed electric aperture field 
and H is the incident tangential magnetic field. To 
extend the formulation over the full range (i.e. to include 
conducting regions), the current densities have to be added 
to equation (2.12) to give: 


e[J(x,y)J 


it inc +^- 

“ JWU 



a 0 
mn mn 

^ mn -k o 



~ a mn 



G 


1 


t 


,exp(j( a m x +e m _y)] 
mn K mn J 


(2.13) 
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A 

where 0 is the complement of the truncation operator 
defined as: 

9[X(?)]=X(r) for r in the aperture (2.14) 

and 9[X(r)]=0 for r in the conducting regions 

and 9 [X( r ) ] =X( r ) -9[X( r ) ] (2.15) 

Note that in equation (2.13) J and E fc are both the unknowns 
to be solved for. 

Equation (2.13) can be recognized as the discrete 

Fourier series for a periodic sequence [25] . Note that 

there is a direct duality between the (x,y,z) domain and the 

(k ,k ,k ) domain. Since all the functions involved here 
x y z 

have a 2tc /m and a 2Tt /ri periodicity in their exponents, one 
period (i.e. one cell) of the structure is sufficient to 
completely specify the transform. That leads to the use of 
the discrete Fourier transform which can be evaluated very 
efficiently by the Fast Fourier transform. It should be 
noted here that because of the exactness of the duality 
between the two domains, no aliasing effects will appear 
when the FFT is performed. By aliasing we mean overlapping 
of spectral components. 

Besides equation (2.13), the boundary condition that 
governs the behavior of the tangential components of the 
electric field, E, over the conducting regions has to be 
satisfied. Equation (2.13) can now be rewritten as : 

= Z _1 F[(-H inC ) + 8(J)] (2.16) 
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where is the Fourier transform of E t 

F is the Fourier transform and F 1 is its inverse 
2, is the product of the Floquet expanding modes 
and Green's dyadic in the spectral domain. 

If the induced currents were available, the solution of "S’ 
could be immediately obtained from (2.16). In practice, 
however, J is an unknown to be solved together with E t and 
equation (2.16) cannot be solved directly. Instead, using 
equation (2.16) a recursive relationship between the (p+l)th 
approximate value of and the pth approximation of E fc is 
now derived and both E fc and J are computed simultaneously, 
via the following iterative procedure: 

a) Start with a guess for E^. in the (x,y) domain and 
apply the truncation operator (i.e. apply the 
boundary condition that E t =0 over any perfectly 
conducting surfaces). 

b) Take the Fourier transform of E^ 

c) Solve for J (p) = F _1 [2F< 9 E fc (p) )]+H t lnC 

(2.17) 


d) Set currents J equal to zero everywhere except over 
the conducting surfaces, that is find: 


A A 

0( j)»e 


F” 1 tz h O E t (p) )]+I t 


me I 


(2.18) 


Substituting equation (2.18) into (2.16) yields: 
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E t (P+ l)=[z"lF 9 F _1 2 F®]E t (p) +Z -1 F [-H t inC +9(H t inC )] 

(2.19 -a) 

And finally taking the inverse Fourier transform of (2.19-a) 
yields: 

E t (p+1 ) =|f~ 1 z" 1 F SF’ 1 ! e E t (p) + F ” 1 ? -1 F[-H t inc+ ®(H t inc )] 

(2.19-b) 


Note that once E fc is evaluated J can also be computed. 
Equation (2.19-b) could be cast in a more convenient form 
(operator form) as: 



E t (p+1) =L l t (p) +? 

(2.20) 

where 

L = F" 1 Z" 1 F^F' 1 Z 

6 is an operator 

and 

C = f 1 *- 1 fl-H t inC 

+9( Hj. ^ nc ) ] is a constant that 


depends on the initial conditions and the incident wave. 


The two most attractive features of this method are the 
following : 

a) No extreme amount of computer memory storage is 
required. 

b) No explicit knowledge of appropriate basis functions 
is needed. 

However, like most iterative techniques, the basic iterative 
scheme suffers from convergence problems. These problems 
and the attempts to alleviate them is the subject of the 
next section. 

2.2 CONVERGENCE OF ITERATIVE SCHEME 

To achieve convergence the important condition that has 
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to be satisfied is that p (L) < 1 or that the spectral 
radius of the operator L has to be less than one. As it 
turns out, for two dimensional cases where the wire spacing 
is greater than two wavelengths, p (L) < 1 and hence 
equation (2.20) converges very quickly for any type of 
incident polarization, angle of incidence and wire 
thickness. However, for spacings less than two wavelengths 
the method fails miserably. To achieve convergence in those 
cases the Successive Relaxation method could be employed to 
"relax" the process and force p (L) < 1 for some relaxation 
factor 0. The choice, not only of the optimum relaxation 
factor, but simply of a relaxation factor that would produce 
a convergent scheme is a difficult task indeed. 

In the one dimensional problem (parallel grid) the 
Contraction Mapping Theory was used very successfully to 
obtain the optimum relaxation factor 0 which forces the 
spectral radius to be less than 1 [26] . To show how this 
theory was used, equation (2.20) is rewritten as: 

g(x n ) = x n+1 = L x n + C (2.21) 

Define a new mapping G(x n ) so that: -• 

G(x n ) =0x + (1-0) g(x n ) (2.22) 

According to the contraction mapping theory [27-31] a 
transformation G of a metric space X onto itself is Lipshitz 
continuous if there exists a p, independent of x and y such 
that 

d(G(x) ,G(y) )<Pd(x,y) for all x,y,e X where d(x,y) is a 
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proper metric in X. For strictly contractive mappings is 
less than one. 

2.2.1 One Dimensional Case 

For the one dimensional case the simplest possible 
metric d that can be used to obtain the optimum e is 
chosen as follows: 

|G(y)-G(x o )| < p |y-x I for P< 1 (2.23) 

Let y=x Q + 6 then 

|g(x q + 6 )-G(x q ) < p |6| or G(x q -*6 )-G(x o ) <p (2.24) 

6 

So the necessary and sufficient condition for contraction 
mapping becomes: 

d (G( x ) ) < p (2.25) 

'dx 

Now substitute (2.22) in (2.24) to obtain: 

|e (x o +6) + (l-6)- g(x o +6)-ex 0 +(l-e) g(x Q )| < P|a | or 

|e+ (1-6) dg(x) I <p 
1 dx ' 

Setting p =0 in the above equation and solving for 6 
yields : 

e = ( dg ( x ) /dx ) / ( dg ( x ) /dx- 1 ) (2.26) 

This value of 0 is called the “contraction" factor since it 
will yield a convergent scheme even in those cases where the 
basic iterative scheme of equation (2.20) fails to converge. 
It should be noted here that in the above analysis e is 
treated as a constant when in fact it is a function of x. 


The reason for that treatment is that 8 is solved in the 
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neighborhood of a solution (root) x Q where the values that 8 
acquires are approximately equal. Therefore 8 can be 
assumed to be constant within that particular neighborhood. 
(For examples and results on the one dimensional problem see 
[26]). 

2.2.2 Two Dimensional Case 

In two dimensions, the basic iterative scheme of 


equation (2.21) is given by: 




Lll x n + L12 ,y n + C1 
L21 x n + L22 y n + C2 


(2.27) 

(2.28) 



four partial derivatives and h y must satisfy the 

following condition [32-35] : 

l 9 xl + |9yl < kl 

| h x|+ |h | < k2 (2.30) 

for kl and k2 less than one and for all points (x,y) in the 


neighborhood R of the root (xo,yo), where R consists of all 
(x,y) with|x-xoj< e , |y-yo|< e , for some positive £. For 
wire spacings less than two wavelengths condition (2.30) is 
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not satisfied. Thus, one has to construct new mappings 
(functions) for the system in (2.27) to obtain convergence 
in a manner similar to the one dimensional case. Now, to 
apply the method of contraction mappings the system (2.29) 
is rewritten as: 

G(x n ,y n ) = 9 x x n + d“© x ) g(x n ,y n ) 

H(x n ,y n ) = 6 y y n + (1-6 ) h(x n ,y n ) (2.31) 

where 6 and 6 are relaxation factors, 
x y 

Unlike the one dimensional problem, Gx, Gy, Hx and Hy 
cannot be separately set equal to zero since they would 
produce a system of equations that are impossible to solve 
for g x =Q, g y =0, h x =0 and h y =0, i.e. 

6 * + (1 -V 0* - 0 

(l-e x ) g » 0 for 8 x (2.32) 

and 

® + (l-e ) h =o 

y y y 

(l-e ) h = 0 for 0 (2.33) 

y x y 

One way to avoid this difficulty is to set kl and k2 to 
nonzero values but their absolute value must always be less 
than one. 
a) First Method 

Let kl and k2 less than one in equations (2.32) and 
(2.33) to obtain: 

I 9 * + (1 -V 9x1 < !/ 2 
|<l-e x > 0 y | < 1/2 


(2.34) 
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and 

| 8 y + d-«y> h y | < 1/2 

1(1-9 ) h I < 1/2 (2.35) 

y x i 

Since hx, hy, gx and gy are complex numbers that implies 
that © x and 0 y can acquire complex values and hence © x and 
0 y are expressed as: 

® x = a + jb ( 2.36-a) 

0 y = c + jd ( 2.36-b) 

Moreover, let 

real (g x ) = a imaginary (g x > = 0 

real (g y ) = y imaginary (g y ) =6 (2.37) 

Upon substituting equations (2.36) and (2.37) into (2.34) 
one obtains: 

|a + jb + (1-a-jb) ( a+ jb) | < 0.5 (2.38-a) 

| (1-a-jb) ( y + j 6) | < 0.5 (2.38-b) 

Taking absolute values yields: 

[(a +a - a a + b0 ) 2 + (b +0 - a& - ba ) 2 ] 1/2 < 0.5 
((Y - 3Y+ b 6 ) 2 + (6 - a 6 - bY ) 2 ] 1/2 < 0.5 
or 

(a +a -a a + b0 ) 2 + (b + 3- a 0 - ba) 2 < ( 0 . 5 ) 2 

( 2.39-a) 

(Y-aY+ b6 ) 2 + ( 6- a 6 - bY ) 2 < ( 0 . 5 ) 2 

( 2. 39-b) 

Equation (2.39) can be expanded to yield two nonlinear 
equations in two unknowns a and b of the form: 
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A1 a 2 + A2 b 2 + A3 a + A4 b + A5 =0.40 2 

a 2 + b 2 -2 a +1 = .4 2 /( Y 2 + 6 2 ) (2.40) 

where Al, A2, A3, A4 and A5 are constants that depend on g x 
and g^. Similarly, to solve for 0^=c+jb another set of 
nonlinear equations is to be solved: 

Bl c 2 + B2 d 2 + B3 c + B4 d + B5 = 0.40 2 

c 2 + d 2 -2 c + 1 . = .4 2 /(e 2 +n 2 ) (2.41) 

where e=real (h ), n=imag inary (h ) and Bl, B2, B3, B4 

X X 

and B5 are constants that depend on h and h . 

x y 

The solution of these nonlinear equations give 9 and 9 

x y 

that would be expected to yield a convergent scheme but, 
unfortunately, they fail to do so for a wire mesh. This 
failure is attributed to the fact that the chosen metric is 
not the appropriate one for this type of geometry, whereas 
it could be a good choice for other geometries of frequency 
selective surfaces. This fact leads to another choice of a 
metric space, 
b) Second Method 

This time the Euclidean norm is chosen as follows: 

| m || 2 = ( |Gx( 2 + |Gy | 2 + |Hx| 2 + | Hy | 2 ) 1/2 

(2.42) 

It is desired to solve for Gx, Gy, Hx and Hy that are 
functions of 9x and 9y with the hope to yield jj M Jj 2 <1* So 
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the basic minimization scheme for solving for x and y in 
this case is the following: 


3 HI 2 

<3© 

w x 


and 


<3 INI 2 
3®y 


0 


(2.43) 


It was found previously that Gx, Gy, Hx and Hy can be 

written in terms of x, y as: 

Gx = 0 + ( 1- 6 )g 

x x' y x 

Gy = U-ejgy 

HX = (1 - e y> h x 

Hy = 0 y + ( 1- ©y ) hy (2.44) 

and hence 


|j = (Gx Gx* + Gy Gy + Hx HX + Hy Hy ) 


(2.45) 

Substituting equations (2.44) into (2.45) and after a long 


and tedious manipulation one obtains: 

2 si. _ m ,Q si. ’ rt /i 1 i 2] 


llM|^=(ai Al+a 1 BHe i Al+j0 1 Cl+ a 2 A2+ a 2 B2+3 2 A2+jB 2 C2+d) 


1/2 


if 

where Al= \g^ +1+ |g x p -g x ~g x 

B1 ’ 9 x - l 9 xf + 9 x' ' I 9 / - 2 I 9 / 


(2.46) * 


Cl=g -g 


A2 *l + |h/-h y - h y *| h / 

B2 *V V - N 2 - 2 N : 


C2=h -h 

y y 


D -IfJ + l h xf + | h yf + l 9 / 


Now to solve for © x = a^+j 3 ^ and 8y = 


a 2 ^2 one 
needs to solve the following system of equations: 
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F 1 (o 1/ e i ,a 2 ,3 2 ) = 


a i ' ^1 ' a 2 ' ^2 ^ = 


M 


d 

■ » 1 

a i 


0 = 

2 

in 2 

[2a 

d 

IM 

2 _ 

n — 

1 

1 

[2 a 2 A 2 +B 2 J 

d 

a 2 


u — 

2 

1 MI 2 

d 

ll»| 

ii_ _ 

n — 

1 

1 

[2 3 1 A 1 +jC 1 ] 

d 



u — 

2 

11 % 

d 

ll M | 

\l_ = 

n — 

1 

1 

[2 3 2 A 2 +jC 2 l 

d 

& 2 


u — 

2 

||M|| 2 


F 3 ( a l r ^l' a 2'^2* = ' 


F4 ( CLj. , 3^ / a 2 ' ^2 ^ ~ 


(2.47) 

By using Newton's method or any other minimization method 

one can solve for , a 2 3^ and 3 2 which will give 

the values for 0 and 0 Unfortunately, once more the 

x y 

values of 0 and 0 obtained by this method yield values 
x y 

| JM j | 2 > 1 for some points inside the cell. It should be noted 
here that the condition that J [ M j j 2 < 1 should be satisfied at 
all sampled points in the cell, and the violation of this 
condition at one point is enough to affect all the other 
points since they are all related together via the two 
dimensional Fourier Transform. 

JPigure (2.3) shows a 16 by 16 array of sampled cell 
points and the value of ||m|| 2 at each point. It can be seen 


that the condition the ||M | 2 <1 is violated at numerous 
points, which implies that a contraction mapping cannot be 
achieved by this method. It was observed that the smaller 

become , 


the wire spacing the larger the values of J | M 
especially near the edges of the wires. 
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Fig . 2.3 The values for || m|L at each sample point inside 
aperture ^ 
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Since neither one of the previous chosen metric spaces 
appear very promising for this particular geometry of 
frequency selective surfaces (i.e. a planar mesh) the trial 
of different metric spaces is put to an end and a different 
line of thought is followed in the next method, 
c) Third Method 

Instead of using 9 and 9 , four different relaxation 

x y 

factors ©11, 912, ©21 and 922 could be utilized to offer 
more degrees of freedom in satisfying condition (2.30). 
Thus, the new modified system of equations becomes: 


x n+1 = 

911 

912 


x" 


( 1-911 ) 

- 912 






+ 



y n+1 = 

921 

922 


_y" 


-921 (1 

- 922) 


g(x n ,y n ) 


=G(x n ,y n ) 


|h(x n ,y n ) 


„ , n n. 
=H ( x ,y ) 


(2.48) 


Now it is easy to set all four partial derivatives Gx, Gy, 
Hx and Hy equal to zero to obtain: 


Gx=911+( 1-911 ) g -912 h =0 

X /v 

Gy=912+( 1-911) g y -912 h y =0 
Hx=921-921 g +(1-922) h =0 

X X 

Hy=922-921 g y +(1-922) h y =0 (2.49) 

Solving this system of equations for 911, 912, 921, 922 
yields : 


911 = 


h g +g ( 1-h ) 
x ^y *x y' 


h g - 
x *y 


(1 " g x ) {1 “V 


(2.50) 
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912 = 


921 = 


922= 


h x g y - d-fl x > (l-h y ) 

h x 

h x g y - (1 ‘ g x ) (1 " h y> 

h g + ( 1-g ) h 
x ^y ^x y 

h g - (1-g ) ( 1-h ) 
x *y y x' ' y ' 


(2.51) 


(2.52) 


(2.53) 


Again, this choice of <9's works very well for the one 
dimensional problem but it does not lead to convergence for 
the two dimensional wire mesh problem. 

To explain why this method does not work for the two 
dimensional problem the theory for constructing convergent 
iterations for a pair of trancendental equations is invoked. 
According to this theory the original system of equation 


n+1 , n 

x =g(x 

n+1 . , n 
y =h( x 




) 

) 


can be written as: 


x n+ ^=x n + a [g ( x n ,y n ) -x n ] + 3 (h(x n ,y n )-y n ]=G(x n ,y n ) 
y n+1 = y n + y tg(x n ,y n )-x n ]+ 6 [h(x n ,y n )-y n ] =H(x n ,y n ) 

(2.54) 


Note the similarity of the above equation with equation 
(2.48). The parameters a, 3 , Y and 6 play the same role in 
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equation (2.54) as the relaxation factors ©11, @12, 021 and 

022 in equation (2.48). To find the root of equation (2.54) 

it is desired to determine a , 0 , y and 6 , by the four 

conditions that the first partial derivatives of G and H are 

zero at some point (x,y) that hopefully is near the root. 

Note that the unknown parameters enter linearly in the same 

way as e's do in equation (2.48), so the calculation of the 

partial derivatives G ,G ,H and H posses no problem. For 

x y x y 

the case of trancendental equations, it is known that this 

method of constructing convergent schemes works provided 

that the partial derivatives g ,g ,h and h DO NOT vary 

x y x y 

very rapidly in the neighborhood of the root ( x 0 *y 0 )* Thus, 

although it is easy to produce a G and an H that are well 

behaved at the root ( x 0 ry 0 ) they may behave quite badly a 

small distance away. If this strategy is to be successful G 

and H must not only have small partial derivatives in some 

region, but this region must also include the desired root. 

For the two dimensional wire mesh it was found that the 

derivatives g ,g ,h and h vary very rapidly, especially at 
x y x y 

points close to the edges of the wire. So this fact, and 
the lack of knowledge of the region within which a root 
exists, causes this method to fail, 
d) Fourth Method 

Finally, another method that could be tried to solve for 
x and y is Newton's method. In this case, we start with the 
basic iterative scheme: 
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x n+1 = 

Lll 

L12 

y n+1 » 

L21 

L22 


n 


n 


Cl 

C2 


which gives: 


n+1 

= Lll 

x n + L12 

y n 

+ Cl 


n+1 

= L21 

x n + L22 

y" 

+ C2 

(2.55) 


Since convergence means that for large n xU — » x n+ ^ equation 
(2.25) can be rewritten as: 


x-Lll x -L12 y - Cl =0 
y-L21 x -L22 y - C2 =0 


(2.56) 


Note that Lll, L12, L21 and L22 are operators so one can 
solve the above equation for a root ( x 0 *y 0 ) by employing 
Newton's method. 

The convergence of this formulation though suffers since 

the derivatives, g ,g ,h and h are much larger than one 

x y x y 

for wire spacings less than two wavelengths or so. This 

fact by itself causes this method to fail. 

Figure (2.4) shows how the relaxation factors 8 and 9 

x y 

contract the basic iterative scheme, but still not enough to 
push the iteration into the region, of convergence. 




Figure 2.4 


Contraction effect of different relaxation 
factors . 



26 


2 . 3 COMMENTS 

It is believed that, unlike the one dimensional problem, 
the two dimensional problem has functions and partial 
derivatives that are very steep, so any method that depends 
in a critical way on magnitudes of derivatives will have 
difficulty to converge. It is also believed that all the 
above mentioned methods for obtaining a convergent scheme 
can be very effectively applied to other geometries of 
frequency selective surfaces, such as an array of metallic 
patches, cross dipoles, circular apertures, etc. 

In conclusion, this method works very well for large 
spacings between adjacent wires without making use of any 
relaxation, contraction or variational factors, but it fails 
miserably to converge for two dimensional problems where the 
mesh spacing is less than two wavelengths or so. 
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3. THE SPECTRAL DOMAIN CONJUGATE GRADIENT APPROACH 

Unlike the previous (S.I.T.) method, in this method the 
induced currents and the aperture fields are solved 
separately. The common features of the S. D.C.G. method and 
the S.I.T. method are that they are both solved in the 
spectral domain and that both make use of the fast Fourier 
transform. In the S. D.C.G. approach, the conjugate gradient 
method is employed to improve upon each previous iterate. 
Hence., the method is basically an iterative technique. 

This part of the the dissertation includes the analysis 
and formulation of the problem in the spectral domain for 
both current densities and electric fields and their 
solution via the conjugate gradient technique. Moreover, a 
number of numerical properties for the conjugate gradient 
method are discussed, and ways of terminating the iterative 
process are suggested. 


3.1 REVIEW OF THE CONJUGATE GRADIENT THEORY 

Suppose that the system that is to be solved is given 

by: 

A "x = y (3.1) 

Let S (0) be and initial guess for x and the residual error 


vector be: 


-*( 0 ) -*■ . -*( 0 ) 
r = y - A x 


(3.2) 


-1 


If A is symmetric positive definite then A is also 
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symmetric positive definite. Now define the quadratic error 
functionals as: 

ERRF X = r A A r =<r ,A 1 r> 

ERRF 2 = r* r = <r,r> =|| r|j 2 (3.3) 

ERRF 3 = r* (A A*) -1 r = <r, (A A*) -1 r> 

where the asterisk * means the conjugate transpose. 

All error functionals in equation (3.3) are positive for all 

values of x^^ except for x^ 0 ^ =x , where jc is the exact 

e e 

solution of Ax=y. In the case where x^^ is equal to the 

exact solution x all the error functionals in (3.3) would 

e 

be equal to zero. 

Now, substitute equation (3.2) in the first error 
functional of equation (3.3) to obtain: 

ERRF 1 =< ( y -Alc (0) ) , A -1 (y - A x (0) ) > or 

ERRF 1 =< ( y , A -1 ~y> - 2<y,x (0) > +<x (0) , a1c (0) > 


ERRF 1 is now a quadratic equation function in x^®^ 


(3.4) 

Let 


:(n) 


be a point in N-dimensional space. Then the equation 


(0) = t (n) + Oj, p 


( n ) 


(3.5) 


■(n) 


is the equation of the line through point x'“' in the 
direction of "p^ n ^ f called the direction vector. For a two 
dimensional interpretation see Figure (3.1). The parameter a 


( 0 ) 


— x 


(n) 


is proportional to the distance jx 
equation (3.5) in equation (3.4) leads to: 


Substituting 
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There are two basic methods that can be used here to 
obtain the next trial vector. The first is the steepest 
descent method and the other one is the conjugate gradient 
method. These methods differ only in the choice of their 
direction vector p^ n ^. Sarkar showed how the steepest 
descent method method can be applied for electrostatic 
problems [36] . In the conjugate gradient method, the 
direction vectors, p^ n ^, must mutually orthogonal with 
respect to the the matrix A. That is, 

<p (l) , A p (:i) > =0 for i H (3.9) 

The iterative scheme of the conjugate gradient method, which 
can now be used to yield successive approximations towards 
the correct solution, is given by Hestenes and Stiefel [37], 
and A, Jennings [38] as: 



(3.10) 
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The equations for the n th iteration are 


s (n) = A £ (n) 


a = 
n 


<p (n) , r (n) > 

<5 (n) , t ln) > 


r(n+l) . ;<n) + ^ *(n> 

:( n+1 ) _ - 2 :(n) _ -|<n) 

n 


(3.11) 


B = 


n 


<; ( n+l> , l(n) > 

<p <n) , 1 (n) > 


t (n+1) = -J(n+1) + -J(n) 


It can be shown [38] that the following othogonal 

relationships are also satisfied: 

-*■( i ) -w i ) 

<r ' ^ , p u '> =0 for i > j (3.12) 


<r 


(i) 


, *<*>> =( 


for i ^ j 


(3.13) 


3.2 CURRENT DENSITY FORMULATION 

The magnetic field H due to an electric current density 
J is given by: 

V x A (x,y,z) 

H ( x ,y ) (3.14) 

U 

-*■ 

where A is the associated magnetic vector potential. A and 
J are related by the free space Green's function 
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A A 

_ exp(-jk.r) 

G(r) 

4 n r 


as follows: 

G ( r , r ' ) . J(r') (3.15) 

From this the electric field intensity E can be derived 
from Maxwell's equations and expressed as: 



E s (x,y,z) = -ju A(x,y,z) + 

joue (3.16) 

For a planar structure we set the z-component of the 
magnetic vector A equal to zero. Now, upon expanding 
equation (3.16) in cartesian- coordinates we obtain, for z=0. 


y^.A(x,y ,z) 


_ E S (x,y) 




/ G * 


Jx 


/g. 


Jy 


(3.17) 


Considering the periodicity of the two dimensional structure 
shown in Figure (3.2) (planar structure), and taking the 
Fourier transform of equation (3.17) leads to: 


E(a ,3 ) = 

mn mn 


jo)U 


. 2 2 
k - a 
o mn 


- a 3 
mn mn 


-a 3 
mn mn 


k 2 - R 2 
k o ^ mn 


— -f 

G J 


(3.18) 

where the sign (") denotes the Fourier transformed quantity. 
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a and 3 represent the Floquet coefficients which were 

defined in the previous chapter as: 

a =2u m/a -ko sindcoscp 
mn 

and 

3 =2ti n/c -2ti m/a cot fl -ko sin © sinm 
mn 


S( a mn' S mn»“ < k ° 2 ‘ “'mn ' 

Fourier transform of Green's function, and J x , 

unknown current densities. 


I is the 
J are the 



Fig. 3.2 Frequency selective surface geometry 
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Notice that the spectrum of 1i s is discrete. That is, it 

exists for discrete values of a and 3 • Note, also, that 

mn mn 

the convolution problem is avoided and instead of dealing 
with an integrodif ferential equation we have to consider 
algebraic equations. 


Taking the inverse Fourier transform of equation (3.18) 
yields : 


2 s <x,y> = 


=-y 

Joe / , 


mn 


k 2 - a 2 - a 3 
o mn mn K mn 


‘Vn 0 


mn 


k 2 - R 2 
k o ° mn 


G . J 


,exp[ + j(a mn x + B mn y)l 


(3.19) 

To enforce the boundary condition over the surface of all 
metallic regions we require that the total tangential 
electric field should satisfy the condition: 

E 1 ( x,y ) + E S ( x,y ) =0 (3.20) 

-►i 

where E is the incident electric field and 

-►S 

E is the scattered electric field 

-V e 

Substituting for the value of E from equation (3.20) into 


equation (3.19) yields 
1 


~E 1 = 


j ooe 



mn 


*o 2 - “ 2 .n 


-a 3 

mn w mn 


-a 3 k 2 - 3 2 

mn M mn o mn 


G ^ a mn '^mn ^ J ^^n '^mn ^ 


.exp (+j ( a mn x+ ^ 


(3.21) 
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Equation (3.21) can be recognized as the inverse 
discrete Fourier transform which can be performed via the 
fast Fourier transform (FFT). Equation (3.21) could now be 
written in an operator form as: 



where Z is the product of G, the Floquet modes and the 
mn 

inverse Fourier transform. 


A solution of the above equation will yield the unknown 

current densities J and J from which the reflected and 

x y 

transmitted fields can be obtained and hence the reflection 

and transmission coefficients could be calculated. It 

should be mentioned here that like the spectral domain 

iteration approach, the spectral domain conjugate gradient 

method is independent of basis functions. 

Now, one way to solve for J and J is to use the 

x y 

conjugate gradient method [19,20,37,38], The conjugate 

gradient method in the spectral domain was used by other 

investigators [21,22] on different geometries. In our case, 

the algorithm of equation (3.11) cannot be directly applied 

on equation (3.22) since Z is symmetric but not self 

mn 

adjoint or positive definite. To over come this difficulty 

and guarantee a convergent scheme, equation (3.22) has to be 

properly modified. To do that, multiply both sides of 

* 

equation (3.22) by Z 

Z ) to obtain: 
mn 


(i.e. the conjugate transpose of 
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-2 E = Z Z J 
mn mn mn 


(3.23) 


where the product Z Z is a Hermitian matrix and 

mn mn 

therefore positive definite. That also means that the 
algorithm (3.11) can now be applied to the transformed 
equation (3.23). In fact, one can apply the previous 
algorithm on equation (3.23) without actually forming 


Z mn Z mn explicitly via the following algorithm [37,38]: 

Let be the initial guess and let the initial residual 

vector be: 

r^°^=Z j ^ ^ + E"*" 
mn 

p<°» = Z* r <0 > - 

mn 


ERRF = 


-( 0 ) 


The equations for the n th iteration are: 


a 


n 



z* ?<"> 
mn 

2 


1 5 <n) 

mn F 

~2 


:(n+l) _ ^-(n) + 


a p 
n 


(n) 


ERRF (n+1) = ERRF (n) 


Z* n r (n) 


Z p 
mn * 


(n) 


(3.24) 


r ( n+ l ) = - a Z p (n) 

n mn e 


& n = 



z* 

mn 

2 


I z* 

1 mn 

2 


p (n+ 1) = z* r (n+1) +& P n) 
mn n v 
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In the above algorithm the root mean square error 

Ik* *ir was chosen as the quadratic functional to be 

minimized. This minimization is also called minimization in 

the range [39] . For a minimization of the functional 
* * —1 1 /2 

r (A A ) r 7 , one could refer to the work done by 

Hestenes and Stiefel, T. Sarkar, J. W. Daniel, T. Cwik and 
Appendix [8.3]. 

3.3 NUMERICAL PROPERTIES OF THE CONJUGATE GRADIENT METHOD 
3.3.1 Singular Operators. 

★ 

Although the transformation Z Z appears to render 

the conjugate gradient method universally applicable for the 

solution of linear operator equations, one must be careful 

of the condition number of Z . If Z is almost singular, 

mn mn 

* 

Z Z will be even more ill-conditioned than Z„ . For 
mn mn mn 

example, let a matrix A be: 



This matrix has a condition number approximately equal to 
400, i.e. 

A.2 / A.]. = 400 

* 

whereas A A has a condition number of = 160,000. That 

means that there is a strong risk of facing poor convergence 
rates. 

One way to check whether or not Z mn is nearly singular 
is to slightly perturb the coefficients of Z mn and apply the 
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conjugate gradient method again. If the results of the 
perturbed system are very different from those obtained from 
the original system, then the matrix Z could be considered 
to be nearly singular and hence poor convergence rates 
should be anticipated. 

3.3.2 Convergence rate 

J. W. Daniel, T. Sarkar and Westlake [40] have shown that 
the convergence rate of the conjugate gradient method is 
given by: 



(3.25) 

where J is the exact solution and X and X are 
e max min 

* 

the maximum and minimum eigenvalues of Z Z in the finite 

mn mn 

dimensional space in which the problem is being solved. In 
this dissertation all problems are solved in a finite 
dimensional space and an investigation of what happens to 
the convergence rate as the dimension of the approximation n 
goes to infinity is avoided. 

W. j. Kammerer and M. Z. Nashed [41] have shown that 
the conjugate gradient method will converge even when Z mn is 
a singular matrix (but with poor rates as mentioned before). 
In that case, the method converges monotonically to a 
solution with minimum norm and the rate of convergence is 
given by Sarkar as: 
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? (n) -J. 


M. 

J (0) + Z + E 1 

mn 


M + n 

J (0) + Z + E 1 

mn 


M 

n 


where M= 


ran 


(2* ) + J (0) + ( Z Z* — ) + I 1 


mn 


ran mn 


(3.26) 


(3.27) 


and Z is the pseudo-inverse of Z 

mn mn 


3.3.3 Stability 

As in most numerical techniques, stability problems may 
appear due to roundoff errors in the calculation of the 
residual and the direction vectors. One possible way of 
automatically detecting instability during the iterative 
process is to look at the ratio a n / a n _i since all scalars, 
cx^, are in the range 

1 < a n < 1 
XTmax X min 

According to T. Sarkar, an upper bound for a n / a n _i is 
X max/ X min and hence stability may be low if X max / X min 
is large. In our case, computational instability with the 
above algorithm for the computation of the residuals was not 
encountered. 

3.3.4 Global versus local convergence 

As with most optimization methods, the conjugate 
gradient method may end up in a local minimum instead of a 
desired global minimum. If the number of unknowns is 
relatively large, it is practically impossible to judge in 
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any way whether or not the minimum found is the desired 
minimum. One possible procedure to use to check this is the 
following: Use several initial guesses in the domain and 

repeat the optimization problem. If all optimizations 
result in approximately the same answer, one could be 
assured that this answer is indeed the desired global 
minimum of the problem. Moreover, it should be mentioned 
here that, from the engineering point of view, one is 
usually not interested in the global minimum if the solution 
obtained can be considered satisfactory. 


3.4 STOPPING PROCEDURES AND INITIAL GUESS 
3.4.1 Stopping procedure 

For the conjugate gradient method there are different 
stopping procedures to terminate the iterative process. 
Three of the most widely used procedures are the following: 


a ) ERROR = 


♦ 

r 


Z J + E J 
mn 


| E 1 


E 1 


< S 


(Normalized error) 


where e is an assigned number of desired accuracy. 

b) Percentage error 

-► 
r 


ERROR% = 


. 100 = 


Z J + E' 
mn 


100 <6 (Normalized error) 


where 6 is another assigned number, 
c) 


2 - 


< C 


ERROR 


(Normalized error) 
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3.4.2 Initial guess 

In all cases checked in this dissertation, the initial 
guess used was = J =0. T ^ at Qi ves a 1/ or 100 

percent, error on the first iteration if the above 
normalized error measures are used. The other reason for 
using a zero guess as a starting point was to see if the 
method converges with the worst possible guess. 

Any other initial guess could be employed to start the 
algorithm. The closer the initial guess to the correct 
answer the better, since the faster the method will 


converge 
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4. FORMULATION OF THE S.D.C.G. METHOD FOR THIN WIRES WITH 
FINITE CONDUCTIVITY AND FOR APERTURE FIELDS 

4.1 EQUIVALENT RADIUS CONCEPT AND INTERNAL INPEDANCE 
The strip analysis can be used to determine the 
scattering characteristics from a mesh of cylindrical wires 
by employing the "equivalent radius" concept. This is 
accomplished by replacing the non-circular cross section of 
a metallic strip with a circular wire whose radius is the 
"equivalent radius" of the non-circular cross section (See 
Figure (4.1)). Butler [42] has shown that the equivalent 
radius of a narrow conducting strip is one fourth of its 
width i.e. 



where a e ^ is the equivalent radius of a cylindrical wire, 
and a is the the width of a thin metallic strip. 



bcO 

r 


Fig. 4.1. 


Equivalent radius of a strip 
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For the case where the wires are of finite conductivity 
the necessary boundary condition that must be satisfied is: 

E 1 + E s = 2. . I (4.1) 

int 

instead of E s + E 1 = 0. where I is the current in the 
wires and Z^ nfc is the internal impedance of the wire which 
is given by Jordan and Balmain [43] as: 


Z int 


m 


2 ua 


eq 


V * a eq> 
r l ( ta eq> 


(4.2) 


where 


m 


/ j (ou ^ //2 
\ m + j^m/ 


is the intrinsic impedance 


1/2 

of the metal. Y is equal to (j u ffl u(o + j.o>e m )) 

and I and I. are the modified Bessel functions which can 
o 1 

be written in terms of infinite series as: 

OO 

- ( x/2 ) 2s+n 

s 1 ( s+n) J (4.3) 

s=o 



n=0 for I and n=l for I, . The case of particular interest 
here, occurs for frequencies sufficiently high that the 
depth of penetration is small compared to the radius of the 
wire. This implies that | Ya e g|>>l and, using the 
asymptotic expansion for I Q and 1^, I o =I^ (See Figure 
(4.2)). Thus, the internal impedance can now be written as: 
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Z. (high frequency) sr 

2 n a eq (4.4) 

2.4 

2.0 

1.6 
1.2 
0.8 
0.4 
0.0 

I 2 3 

Fig. 4.2. Modified Bessel functions 

From equation (4.4) it can be seen that for a small skin 
depth Z^ nt is equal to the surface impedance of a thick 
metal sheet that is one meter long and 2n a g meters wide. 
Now substituting the expression for in equation (4.4) 



yields : 
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This expression for Z^ nfc can now be used in equation (4.1)/ 
i.e. 


1 - (Z lnf A) J 


since J=I/A where A is the surface area of the wire. 
This leads to: 



Now equation (4.8) can be solved for J using the algorithm 

mentioned before in equation (3.24) and replacing Z mn by 

(Z — Z • . ) . Rather than form the matrix ( Z -Z . . ) 

mn mt mn mt 

explicitly, one can carry out the calculation using the 
following algorithm: 



46 



The equations for the n th iteration are: 



“(n+1 ) _ * -►( n+1 ) =* ♦ ( n+1 ) + B $(n) 

p ~ Z ran r Z int r n 

END OF DO LOOP 


4.2 SOLUTION OF APERTURE FIELDS 

To solve for the aperture fields (See Figure (4.3)), 
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equation (2.12) is used as the starting point, i.e. 


H 1 =• 


2j 


UU 



mn 


a 3 ko 2 - a 2 

mn p mn mn 


~ ^ mn'^o ^ ~ a mn ^ mn 


G( a ,3 ) E a 

mn mn 


* expC+j(a mn x+p mn y )] 


(4.10) 


Fig. 4.3. 


Sampling for the Aperture fields 



For a complete derivation of equation (4.10), see Chapter 2 

and Appendix 8.2. H 1 is the incident magnetic field which 

is a known quantity, G and a ,3 were defined before 

mn mn 

and they are also known. The unknown in this case is E a , so 
equation (4.10) can now be cast into operator form as: 
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= Y mn 13 (4.11) 

Y mn ' like Z mn in Chapter 3, is neither positive definite nor 

self adjoint/ so both sides of equation (4.11) are 

multiplied by the conjugate transpose of Y , i.e. 

' mn 


Y H = Y Y E 

mn mn mn 


(4.12) 


Now the algorithm in equation (3.24) can be applied to the 
above equation be replacing Z mn by Y mn / J by 1 a , and -S 1 by 


As before with the current densities, we choose as 

an initial guess for the aperture field, E a , and start 
iterating. In this dissertation the initial guess is chosen 
to be equal to zero in all check cases and this leads to the 
following stopping procedure: 


ERROR 




(4.13) 


If a percentage error is desired then the stopping procedure 
becomes : 


ERROR% = 


H 1 - 


■iXlOO = - 


' 

H 1 


mn 


x 100 


(4.14) 


Note that, for a zero initial guess, the first error 
estimate will be equal to 1 (for the first iteration), 
whereas, the second estimate will yield a 100% error. 



49 


4.3 REFLECTION COEFFICIENTS 

The transmission and reflection coefficients are the 
quantities of most important in characterizing the 
properties of a mesh. In order to define those coefficients 
for both polarizations, transverse electric (TE) and 
transverse magnetic (TM), it is necessary to first define 
the incident and scattered fields. 

For TE. polarization , the incident fields are: 

E x = E q sin(-cp) ; Ey = E Q coscp 


E coscp cosd 

Hx . • 


E q sinco cosd 


where E q is the amplitude of the incident electric field and 
1 /2 

h= ( uo/eo ) ' is the free space wave impedance. 

For TM polarization, the incident fields are given by: 

E x = E q cos& coscp ; E y = E q cosO sincp 


H = 
x 


E q sin ( cp -ix /2) E Q cos(cp -n /2) 

■ ■ — ; H = 

y 


h J n 

According to Wait and Hill [44] , when the spacing 
between adjacent wires of the mesh is less than X/2, there 
is only one grating lobe and only the constant current 
components J QOX and J Qoy contribute to the scattered field. 

J oox an< ^ J ooy are zero " mo ^ e current density components. 

The rectangular components of the scattered field, for large 
z , are given by: 



t tu 
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s 2 2 2 

E x - J oox^ 1 *“ s ^ n ^ cos v ^ -J ooy si,n ^ s i n <p coscp 

exp ik [zcosd + sind (x coscp + y sirup)] 

(4.17) 

s 2 2 2 

E * J qox s * n ^ s i n «P coscp “J 00 y (l.-sin d sin cp) 

exp ik [zcosd + sind(x coscp + y sincp)] 

(4.18) 

The above expressions can also be obtained from equation 

(3.18) as follows: Solve for and substitute the solution 

in equation (3.18) to obtain the scattered fields. That is. 



and so the reflection (amplitude) coefficient becomes: 


R = E S / (E 2 + E 2 ) 1/ ' 2 

x x ' x y 

R = E S / (E 2 + E 2 ) ly/2 

y y x y 


(4.20) 

(4.21) 


If. the total power reflection coefficient R is desired 
then the following expression can be used: 
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where E is the scattered field due to the induced current 
densities, J, derived from equation (4.19) and, after taking 
the inverse Fourier transform, H S is the scattered magnetic 
field derived from by making use of Maxwell's equations. 

Moreover, if the total power transmission coefficient, 
T, is to be computed/ one can employ the formula below: 


Real 


x H 


" a * (- z) dA 


(4.23) 


aperture- 


Real 


aperture 


/ 


E x H .(- z) dA 


where E a is the aperture Electric field and 15 a is the 
magnetic field in the aperture derived from Maxwell's 
equation 


-j o)U H a =V x E a (4.24) 
For perfectly conducting frequency selective structures it 
is also true that: 



this condition can be used in the perfectly conducting cases 
to check the convergence and accuracy of the results. 
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5. RESULTS AND COMMENTS 

5.1 ONE DIMENSIONAL CASE (INFINITE GRATING OF PARALLEL 
STRIPS) 

5.1.1 Current densities 

The one dimensional case was studied first, and 
compared to the Spectral Domain Approach with the 
contraction factor [26] , since results from that method were 
readily available. In Table (5.1), for example, the current 
density levels, for thin strips obtained by both the S.I.T. 
method and the S.D.C.G. method are in very good agreement. 
This implies that either method can be employed to generate 
the induced current densities on a strip or wire grating for 
any incident field. 


Table 5.1. Current densities obtained by the S.D.C.G. 

method and the S.I.T. method. (See Figure (2.1) 
for the geometry), a is the spacing between 



adjacent 
strips . 

strips and w 

is the width 

of the 

Spacing (a) 

Width(w) 

S.D.C.G. 

S.I.T. 

Difference 

0.55 A 

0.005 A 

0.02664928 

0.02770429 

0.001055 

0.25 X 

0.005 A 

0.05155611 

0.05183827 

0.000281 

0.125 A 

0.005 A 

0.07172995 

0.07114100 

0.000588 

0.100 A 

0.002 A 

0.07545375 

0.07521373 

0.000240 
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Figures (5.1) and (5.2) show that the current densities 
obtained via the S.D.C.G. method for a parallel grid with 
thick strip and a normally incident field behave as 
expected. For the copolar component the current density 
curves downwards. That is, it is larger at the edges than 



Fig. 5.1. Amplitude of y-component of the current density 
for a grid of parallel strips and for a normally 
incident field ( 0=0°). The incident electric 
field is along the y axis. 
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at the center (Fig. 5.1), a phenomenon attributed to the 
edge effects of the metal strip. On the other hand, the 
cross-polar component in Figure (5.2) curves the other way 
around, i.e. outwards. 



Fig. 5.2. Amplitude for x-component of the current density 
for a grating of parallel strips and for a 
normally incident field ( S = 0°). The 
polarization is TE and the incident electric 
field is y directed. 
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Figure (5.3) shows how the reflection coefficient for normal 
incidence increases as the spacing of the grid gets smaller. 
This anticipated behavior is due to the fact that the closer 
the wires, the closer the grid structure resembles a solid 
metal sheet. Notice that this method is even valid for 
spacings of 1/100 of a wavelength; a fact that renders this 
algorithm very useful for radiometric applications where the 
spacing between wires is of the order of 1/10 of a 
wavelength or less. It should be mentioned here that the 
data in Figure (5.3) are not compared with any measured data 
or calculations made using other methods since at these 
spacings neither calculations nor measured data exist. 



Fig. 5.3. Reflection coefficient for a grid of parallel 
strips and spacings of 1/10 X and less. 
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The only other method that can generate reflection 
coefficients at those spacings is the S.I.T. modified with 
the contraction factor given by [26]. Table 5.2 shows that 
the reflection coefficients obtained by both, the S.D.C.G. 
method and Brand's method [26] are almost identical for 
various wire spacings. 


Table 

5.2. 

Reflection coefficients for 
spacings calculated by the 
the contraction factor-S.I. 
incidence ) 

different wire 

S. D.C.G. method and 

T. method. (Normal 

Spacing 

S.D.C.G. 

S.I.T. 

0.125 

A 

0.844 

0.843 

0.10 

A 

0.888 

0.885 

0.06 

A 

0.954 

0.960 

0.05 

A 

0.967 

0.969 

0.02 

A 

0.994 

0.994 

0.01 

A 

0.999 

0.999 


The S.D.C.G. algorithm for one dimensional cases (i.e. 
parallel wires) converges in at most six iterations with a 
normalized error of less than 0.5 percent. The CPU time 
used for each of the above cases was about 20 sec on the 
3081 IBM system and for a 32x32 sampling rate. This time 
includes plotting time. Table 5.3 shows how the normalized 
error decreases at each iteration for spacings of 1/10 of a 
wavelength or less. 
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Table 5.3. Normalized error versus number of iterations for 
spacings less than 1/10 of a wavelength 

Spacings between strips Normalized percentage Number of 
(width of strips=0.002X) error ((| r||/||E l ||)x!00 Iterations 

100 1 

0.10 x 13.57 2 

0.348 3 

100 1 

0.09 X 17 2 

0.14 3 

100 1 

0.07 X 29 2 

0.3 3 

100 1 

0.05 X 52 2 

0.3 3 

100 1 

0.04 X 69 2 

0.13 3 

100 1 

0.03 X 69 2 

0.25 3 

100 1 

73 2 

0.01 X 40 3 

17 4 


0.5 


5 
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between the two methods is very good indeed. To actually 
see how close the numbers are, Table 5.4 gives the values 
for the aperture field at each sampling point for both 
methods. 

Table 5.4. Values of aperture field at each sampling point 
for an aperture size of 0.25^. an d normal 
incidence . 


Cell point 

S.I.T. 

S.D.C.G. 

on x-axis 

method 

method 

-0 . 12 9126 31 0 

0. 182 2 59561 E-0 ' 

0. 0000 0 OOOOE* 0 0 

-0 . 120795548 

0. 3840 13295 

0.37743 0677 

-0.112464-»€6 

0.560089946 

0.556 849957 

-0 . 104174083 

0.661797166 

0.660453684 

-0 *95 807 320 4E-0 1 

-0.7 38 4 769 92 

0. 7385-: 7 141 

- 0 • 9 7 47 26176E-01 

0 .796682556 

0 .7978 7993 4 

-Q.791419148E-01 

0. €44159245 

0 .8462 64482 

-3 . 70 til 1S2SE-Q1 

0 .88262^010 

0.835471702 

-0.624S04755E-01 

0. 514659 142 

0.518128343 

-0.S414^746?E-01 

0.940926552 

0.944895143 

-0. 458190 IgO E-01 

0. 5625 533 82 

0.966925740 

“0 • 274 9 92 85 ? E— 0 1 

0.579851842 

0.9R454 79 73 

-0 .29157558 6 E-0 1 

0.593371725 . 

0. 998 3290 93 

-0.208 26 8240E-01 

1 .00 326443 

1 .00839133 

-0. 12496091 4E-0 1 

1 .00979042 

1 .0150 3463 

-0.416536257E-02 

1 .01 3 01479 

1.01830196 

0 .416536257 E— 02 

1 . 01 301575 

1 .01 8310 55 

3. 124 96 0 91 4 £-.11 

1 .009 76042 

1 .01503161 

0.20924 e240E-01 

1 .00 326347 

1.00839233 

0.291575586E-01 

0. 593371725 

0.998323321 

0.374 8 82 8 5 5 E-0 1 

0 • 5798 513 42 

0 .984559417 

0 .458190 180 8-01 

0 .562S536 55 

0.56691 4892 

0 .5414«»74€ 8 E-01 

0.940926552 

0. 94490 2539 

0 .624604795E-01 

0. 514659023 

0 .918127894 

0 .70611 152SE-01 

0.882 627010 

0.8854 74324 

0. 791419148E-01 

0. €441 59565 

0. €46 2 53 52 1 

0.87 472 6 176E-01 

0.79663 2477 

0 .7978348 22 

0 .95 € 03 22 0 «£- 0 3 

0. 7 38 4 76 753 

0 . 738585949 

0.104 17 4 3 P 2 

0 .661 7 >76 43 

0.660460711 

0.112464786 

0.560 0907-0 

0.5568441 15 

0 . 120 7<=55 46 

0.2 €40 14010 

0.377436161 

0.129126310 

0. 1 82273 380E-0 1 

o.ooooooooo=>oo 


Figure (5.5) shows how the error is reduced at each 


iteration for this case. It should be mentioned here that 
this normalized error (See Section 3.4 for definition) is 


the total error one obtains by sampling the unit cell by a 
rate of 32x32 samples. Another check was obtained against 
published data given by Mittra and Tsao [16] . Again, the 
agreement between the two techniques is shown in Figure 
(5.6) . 


NORMALIZED ERROR APERTURE FIELDS ( TE ) 


nuwer of iterations 


The normalized error for an aperture field of 
0 . 25 A. in size 
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Fig. 5.6. Amplitude of the aperture electric field for a 
unit cell with a=1.4A and b=0.6 a . The 
incident field is at normal incidence and with 
TE polarization 


For the infinite grid of parallel wires, the current 
densities, the aperture fields and the reflection 
coefficients obtained by this algorithm are in very good 
agreement with the S. I .T. -Contraction method and the results 
published by Mittra and Tsao. 
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5.2 TWO DIMENSIONAL CASE (i.e. INFINITE GRID WITH SQUARE 

OPENINGS) 

5.2.1 Current densities 

For the two dimensional case a number of cases are 
checked against calculations by Wait and Hill [ 6 , 44 , 45 ] and 
Kontorovich, Astrakham and their colleagues [9,10]. 

Overall, very good agreement is found in the calculation of 
the reflection coefficients for different angles of 
incidence, polarization and wire spacing. The reason for 
comparing reflection coefficients with other methods is 
simply that the reflection coefficient is the parameter of 
most importance in designing wire meshes. 

Figure (5.7), (5.8), (5.9) and (5.10) show these 
reflection coefficients for both transverse electric (TE) 
and transverse magnetic (TM) incidence. Calculations using 
the S.D.C.G. method are compared with two other methods. 
Wait's method is based on a Fourier series expansion 
solution, whereas, the Kontorovich-Astrakham method is based 
on the averaged boundary condition technique. In all those 
figures, a=b represents the wire spacing of the square mesh 
(See Figure 2.1) and c is the equivalent radius of the 
strips. Figure (5.7) exhibits the characteristic Brewster- 
angle minimum for the S.D.C.G. method and Wait's method. 

The discrepancy between the two curves is attributed to the 
fact that in the S.D.C.G. method planar strips are actually 
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used instead of round wires. The sampling rate used in 
these cases was 16x16 samples. For thin wires, this 
sampling rate is good enough to obtain a good estimate for 
the reflection coefficients; this is evident from these 
figures. If more accuracy is desired, the number of samples 
can be increased. In Figure (5.10) one can see that, by 
increasing the sampling rate, a slightly better estimate can 
be obtained. 



Fig. 5.7. Reflection coefficient for TM incidence and 
various spacings for d=70 deg., and cp=0 deg. 




reflection coefficient 


(j )=0° 

3 - = 0.25 2~ b 


x---x * KONTOROVICH 

astrakham. 


02 / 


90 so" 


5.9. The reflection coefficient for a spacing of 

a=0.25Aand for different values of the angle of 
incidence theta. (TM polarization) 


REFLECTION COEFFICIENT 
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Fig. 5.10. The reflection coefficient for TE polarization 

and different angles of theta. The wire spacing 
is equal to 0.25X . 



6 : 


Figures (5.11) and (5.12) confirm the expected result 
that the wider the wires the larger the reflection 
coefficients. This result should be anticipated since a 
mesh with wide wires is a better approximation to a 
continuous metal sheet. 
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Fig. 5.11. Reflection coefficients for different 

thicknesses and for different angles of 
incidence theta. The polarization is transverse 
magnetic and the mesh opening is a=b=0.25X . 
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Fig. 5.12. Reflection coefficients for different widths and 
for different angles of incidence theta. (TE 
polarization) . 


Figures (5.13) and (5.14) depict the change in the 
reflection coefficient when the wire mesh consists of wires 
with finite conductivity. The figures confirm the fact that 
the reflection coefficient of a lossy wire-mesh is less than 
that of the perfectly conducting wires case. The reason for 
this difference is that, for perfectly conducting wires 
( o=oo), the reflection and transmission coefficients are 
governed by the relation: 
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M 2 +m 2 =1 

where T is the transmission coefficient and R is the 

reflection coefficient. For lossy wires, due to the loss of 

2 2 

energy in the wires, T + R #1. 



er=5.l0 2 7 J m 


-i 
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Fig. 5.14. The reflection coefficient for different 
conductivities . 

So far, we have only discussed the reflection coefficients 

for different polarizations, angles of incidence, widths and 

wire spacings. We have also compared them with whatever 

data were available to us [6,9,10], In the figures to 

follow, the current densities are presented and analyzed for 

different cases of interest. First, we start with Figures 

(5.15) and (5.16) where the current densities, J and J , 

x y 

are depicted. The spacing used in that case was 0.25X 
wavelengths and a thickness of 0.005 X . The sampling rate 
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was 32x32 and the wave was normally incident for a TE 
polarization. For wider strips Figures (5.17) and (5.18) 
show how the current densities behave. And for a case with 
lossy wires Figures (5.19) and (5.20) give the results. In 
those cases , the square-shaped unit cell was used. 



Fig. 5.15. Amplitude of the x-component of the current 

density for a normally incident wave on a square 
mesh with a spacing of 0 . 25 X between adjacent 
wires. (Incident E field is along the y axis). 
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Fig. 5.16. Amplitude of the y component of the current 

density for a normally incident wave on a square 
mesh with a spacing of 0.25 X. (E incident is 
along the y axis ) . 
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32X32 TE,TH==70,PHI = 0 (AA=.355 BB = . 25.) 



Pig. 5.17. Amplitude of x component of the current density 
for an obliquely incident wave 0=70°) on a 
square mesh with wide metal strips 
(width= .105 X). TE polarization with E 
incident along the y axis). 
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32X32 TE,TH=70,PHI=0 (AA=.355^ BB = .25^) 



Fig. 5.18. Amplitude of the y component of the current 

density for a wave incident at an angle theta=70 
on a square mesh of wide metal strips 
(width=0.105 A), (E incident is along the y 
axis) . 
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SI.GMA-500, TE.TH=0,PHI=0 (A=.25 


0.005 THICK) 


0.0927 


0.0018 


0 0009 


0 0099 


0 i3 . 



Fig. 5.19. Amplitude of x component of the current density 
for a square grid of thin strips but with a 
conductivity of a- 500 y/m . A normally incident 
field ( 9=0 0 ) and a sampling rate of 16x16 
samples are used. 


76 


SIGMA=500, 


TE.TH=0, PHI-0 (A 25 


0.005 THICK) 



Fig 


5.20 


f the current 

Amplitude of the y|°®P°^ e ^ lth wide metal^^ ^ 

density { The inductivity, is ^/incident 

along the y axis. 
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It should be mentioned here that the magnitude of the 
current densities becomes smaller as the conductivity of the 
metal strips or-wires is reduced. This result should be 
anticipated since the smaller the conductivity of the wires 
the lossier they are. 

Now, to illustrate the significant effects that occur 
at a bonded junction, the cross-shaped unit cell is used. 

The current densities obtained in this case are depicted in 
Figures (5.21) to (5.26) for different spacings, widths and 
angles of incidence. It can be seen from all figures that 
this method predicts the step discontinuity at the bonded 
junction. It should be stressed here that in this 
dissertation only the bonded case is treated; that is, the 
case where a perfect contact between the wires exists at 
each junction. 

Since the existing mesh surfaces resemble more closely 
the bonded mesh, than the unbonded case, a study of the 
unbonded mesh was not done here. Quite often, in practice, 
the wires are soldered at the bonds to obtain a perfect 
contact. The study of the unbonded mesh is of interest, 
though, because of its physical analogy with a thin 
magnetized plasma. Anisotropic unbonded wire mesh can be 
used to simulate a thin magnetized plasma sheet. Wait has 
calculated the reflection and transmission coefficients for 
the unbonded mesh case [6] . 
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Fig. 5.21. Amplitude of y current density component for a 

normally incident wave with TE polarization on a 
square mesh (a=b=1.25X and width =0.02 \ ). (E 
incident is along the y axis). 
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TE,TH=00;PHI=0 (A=B = 1.25^AND .02jj THICK) 


0.63 


Fig. 5.22. 


Amplitude of x current density component for 
normally incident wave on a square mesh. 
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TE,TH=0,PHI=0 (A=. 25^ 0.005j|THICKj 32X32 



0-13 


Fig. 5.23. Amplitude of y current density component for 
normally incident wave on a square mesh. ©=0 
and E incident is along the x axis. 


O CD 
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TE,TH=0,PHI=0 (A=.25^ 0.005^ THICK) 32X32 



Fig. 5.24. Amplitude of x current density component (cross- 
polar) for a normally incident wave on a square 
mesh with thin strips. 
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TE,TH=60,PHI=0 (A=B= 25^AND 


.02^ THICK) 



Fig. 5.25. Amplitude of y current density component for an 
obliquely incident wave on a square mesh with 
strips of width equal to .02 \ . 
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TE,TH=60,PHI=0 (A=B = .25j^AND .02jTHICK) 



Fig. 5.26. Amplitude of x current density component (cross 
polar) for a square mesh with an obliquely 
incident wave. 


Figures (5.27) to (5.30) show how the normalized error 
is reduced at each iteration. It can be seen from all these 
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figures that the residual error decreases monotonically . 

From Figure (5.28)/ one can see that the closer the strip 
spacing, the longer it takes to converge to a specified 
normalized error. The difference in the normalized error 
between the 0.25A and 0.75A spacings is indeed large, 
whereas the corresponding difference in the normalized error 
between the spacings of 0.75X and 1.25A is not that drastic. 



Fig. 5.27. Normalized error for currents for a square mesh 
with a=b=1.25A , theta=30° and phi=0°. 
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Fig. 5.29 
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To see how this occurs, we recall the expressions for 

= "j 


a , B and G = — which are the elements 

mn mn 2 - y / k ^?^ 2 

that form the entries of matrix Z . These elements are 

mn 

functions of angles theta (0) and phi ( cp ) and of the wire 
spacing. This means that any change in theta, phi or the 
spacing will yield a change in the matrix Z mn , and hence, 
the eigenvalues of the matrix will be different. It was 
mentioned before, in Chapter 3, that the rate of convergence 
depends on the eigenvalues of matrix Z mn . Therefore, any 
change in theta, phi or in wire spacing will change the rate 
of convergence. 

Another interesting phenomenon is observed in Figure 

(5.30) where the normalized error for the same wire spacing 
and the same incident field, but for differently shaped unit 
cells is plotted. From that figure it is clear that the 
normalized error for the cross-shaped unit cell decreases 
much faster than that of the square-shaped unit cell. 
Although both unit cells generate the same currents and 
reflection coefficients, the cross-shaped unit cell can be 
more advantageous as far as computing time is concerned. 

One reason for this difference between the two unit cells 
lies in the fact- that in the cross-shaped unit cell the wire 
strips appear to be wider to the algorithm than the 
corresponding strips in the other unit cell, as Figure 

(5.31) illustrates. 



»o»3om omNHrs2»Q2 
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Fig. 5.30. Normalized error for two differently shaped unit 
cells . 
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Before we present some results for the aperture fields, 
the problem of evaluating the reflection coefficients from a 
frequency selective planar surface, shown in Figure (5.32), 
is discussed. Table (5,5) gives the results for the 
reflection coefficient evaluated by this algorithm for 
different values of Q. Here, it is observed that a cross- 
polar component arises even for a normally incident wave. 
This result is very important in assessing the degree of 
depolarization from such a planar structure. 



Fig. 5.32. A different frequency selective surface 
geometry. 

It should be mentioned here that this configuration 
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offers a better approximation to the knitted mesh than the 
infinite square grid. The reason for this is that the 
periodicity of the above planar structure resembles that of 
the knitted mesh. 

Table 5.5. Reflection coefficients from the arrangement in 
Figure (5.32). (Normal incidence and TM 
polarization) . 


» 

copolar 

cross-po 

90° 
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0.0 

00 
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0.16 

60° 

0.7316 
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o 
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5.2.2. 

Aperture Fields 



Figure (5.33) shows how the aperture field is compared 
with the results published by Tsao and Mittra [16] . This 
happens to be the only available data for aperture fields 
that we can compare with our calculations. For spacings 
larger than one wavelength there are more than one 
propagating modes (i.e. whenever k Q > a mn + which 

appear as lobes in the aperture field. Notice the four 
lobes in Figure (5.33) for a spacing of four wavelengths 
between the adjacent strips. Figures (5.34) and (5.35) 
depict the amplitudes of the x- and y- components of the 


c.- 
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aperture electric field for a normally incident field on a 
square mesh with the dimensions a=b=1.25X . Note again, 
that this algorithm can predict the two propagating lobes 
and the edge effects on the strips that are perpendicular to 
the y-directed incident electric field. Moreover, Figure 
(5.36) and (5.37) give the amplitudes of the aperture fields 
for a different type of polarization (Transverse magnetic or 
TM) and a mesh with innerspacing given by a=b=0.25X . In 
this case the angles theta and phi are both equal to 30 
degrees. In Figures (5.38) and (5.39), a smaller spacing is 
used (a=b=0.125 ) and an angle of incidence equal to 30°, 

to see if the algorithm can still converge under the 
conditions of oblique incidence and smaller spacings. 
Finally, in Figures (5.40) and (5.41) a sampling rate of 
16x16 samples is used instead of 32x32. The x and y 
components of the electric field in the aperture are shown. 
In this case the wave is normally incident on a mesh of thin 
strips and s spacing equal to 0.25 X. Once more, the edge 
effects become very evident. 
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APERTURE FIELD TH=PH=0 TE (A=1.25^) 



Amplitude of y component of the aperture field 
for a normally incident wave on a square mesh 
with a=b=1.25A. (Two lobes). 
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APERTURE FIELD TH=PH = 0 TE (A=1.25j>) 



Fig. 5.35. Amplitude of x component of the aperture field 
for a normally incident wave with TE 
polarization . 
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* 


Fig . 


APERTURE FIELD TH=PH=30 TM (A=0 25j0 



5.36. Amplitude of x component of the aperture field 

for a wave incident on the square mesh at angles 
theta=30°, phi=30° and with a TM polarization. 
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APERTURE FIELD TH=PH=30 TM (A=0 25j) 



Fig. 5.37. Amplitude of y component of the aperture field 
for an obliquely incident wave on a square mesh 
with a=b=0.25X . 
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APERTURE FIELD TH=30 PH=0 TM (A=0.125j) 



Fig. 5.38. Amplitude of x component of the aperture field 
for a wave incident on a thin strip mesh with 
a=b=0 . 125 A. . 
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APERTURE FIELD TH=30 PH=0 TM (A=0.i25^) 



Fig. 5.39. Amplitude of y component of the aperture field 
for a wave incident at angle theta=30° on a 
square mesh with a=b=0.125\. 
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TE, NORMAL (A=B = .25^AND 005jJHICK) 



Fig. 5.40. Amplitude of y component of the electric 

aperture field for a normally incident wave on a 
mesh and sampling rate of 16x16 samples. 


TE.NORMAL (A=B= 25 AND .005 THICK) 



Fig. 5.41. Amplitude of x component of the electric 

aperture field for a sampling rate of 16x16 
samples and a normally incident wave. 


102 


5.3 CPU TIME AND STORAGE REQCJIRMENTS 

In general the two dimensional problem takes longer to 
converge than the one dimensional case. One of the main 
reasons for that is the fact -that a two dimensional FFT is 
used and more sampling points are required in the two 
dimensional case. Moreover, in the two dimensional case, 
one has to solve for far more unknowns than in the one 
dimensional problem. This number of unknowns also affects 
the storage requirements for the two dimensional problem. 

In general the CPU time depends on the sampling rate more 
than on anything else. For example, for a sampling rate of 
16x16 samples, it takes anywhere from 30 seconds to 1.40 
minutes of CPU time (on an IBM 3081 system) to converge to a 
reasonably accurate result. This time includes sorting and 
plotting of data. For a 32x32 sampling rate the CPU time 
is, as expected, much more. In fact, in this case the range 
is somewhere between 1:32 and 8:00 minutes. As mentioned 
before, the CPU time also depends on the angle of incidence 
and the strip spacing. 

The program size is 56,064 bytes for a 16x16 sampling 
rate and 149,192 bytes for a 32x32 sampling rate. 
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6. COMMENTS AND SUGGESTIONS FOR FUTURE RESEARCH 

Here, a number of recommendations for future research 
related to the mesh problem and the S.D.C.G. method are 
mentioned. 

One). Skew-Symmetric Configuration for a mesh 




(a) 


(b) 


Fig. 6.1. Different sampling patterns (a) rectangular 
(b) Skew- symmetric . 

For a rectangular or a square grid, the number of 
samples, which is also the number of Floquet modes, 
corresponds to the number of couplings being taken into 
account. Moreover, for rectangular sampling, the Fast 
Fourier Transform can be used directly. On the other hand. 
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for non-rectangular sampling, such as that shown in Figure 
6.1 (b), the number of samples may not correspond exactly to 
the number of couplings being taken into account. That will 
yield some erroneous results. Moreover, in this case FFT 
can not be used and hence a Discrete Fourier Transform has 
to be employed instead. So the actual problem here is to 
modify the existing algorithm in order to represent the 
strips and their width as accurately as possible. The 
reason for doing this is to avoid any aliasing problems that 
may arise from sampling such a configuration. Solving this 
problem is important because a study of the reflection 
coefficients as a function of the angle ¥ will give new 
insight in designing mesh structure that are skew-symmetric. 
Moreover, this configuration might offer a better 
approximation to the actual woven structure than the 
rectangular mesh. 

Two) . Double screen 

The scattering properties of such a structure are of 
interest because a double screen can be used as a filter in 
microwave applications. To solve this problem, the original 
structure is divided into two substructures, and the 
principles of equivalence and superposition are used to 
obtain the formulation in the spectral domain. According to 
Tsao and Mittra the problem of the double screen can be 
represented as in Figure (6.2): 
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Conductor Conductor 

Fig. 6.2. Equivalent problem for a double screen 


The equations for the aperture field, for example, are 
similar to those used in our work for the single square 
mesh, but with the phase difference between the two meshes 
taken into consideration. Tsao and Mittra give the 
following equations in the spectral domain for E a ^ and E a 2 : 

5.E a 1 exp [ j ^ a mn x+ ^mn^ ^ 

1 -*• 

= — H 1 for the electric conductor 



2 
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„ 2_ 2 

o a mn 
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for the magnetic conductor 
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Y mnd 
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Y mn - 


‘I V k o ! ‘ a 2 -n' B 


mn 


2 + 8 2 -k 2 

mn mn o 


2 2 2 
for k z > a + 8 mn 
o mn mn 


for k 2 < a 2 + 8 2 mn 

o mn w mn 


The total aperture field at any point is the superposition 

4a *4a 

of E x and E 2 . 

Three). A mesh over a ground plane. 



MESH 


GROUND PLANE 


Fig. 6.3. A mesh over a ground plane 
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In this problem, the evanenscent field from the mesh 
interacts with the adjacent ground plane, and hence, the 
total scattered field is different than that of the mesh in 
free space. To solve this problem the superposition 
principle should be used as follows. First, solve for the 
field reflected by the ground plane (or dielectric plane). 
Second, find the scattered field due to the mesh and finally 
find the total scattered field. Once the total scattered 
fields are known the reflection coefficients can be 
determined . 

Four). Determination of the Electromagnetic properties of a 
mesh with wires made of different alloys for 
Radiometric operation 

Figure (6.4) shows that the wires used in constructing 
the actual mesh are made of molybdenum 1.2 mill in diameter 
and there is 4-6% gold, by weight, plated over the 
molybdenum. In thickness this corresponds to 8-11 micro- 
inches . 



Fig. 6.4. Gold plated wire substrate 
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The resistivity ( p ) of gold and molybdenum is 

different. For gold p is 2.35 -cm, whereas for 

molybdenum the resistivity is 5.2 uQ~cm. In our work we 

assumed the mesh was constructed with wires of one type of 

resistivity and not alloys. If the skin depth, for a 

certain frequency, is larger than the thickness of the fold 

that covers the molybdenum wire, the field would penetrate 

into the molybdenum region. This means that the wire cannot 

be considered as uniform any more. Therefore, the current 

algorithm has to be modified to take into consideration this 

difference in resistivity. One way to do that would be to 

derive an impedance expression for thin wires made of any 

kind of alloys. Once this impedance is obtained, it can be 

used in the S.D.C.G. method as follows: Start with the 

equation for the currents in chapter 3. That is, 

Z mri J = E S (6.1) 

mn 

The new boundary condition becomes: 

r • A alloy 7 t6 ‘ 2) 


where z a ^i 0 y is the internal impedance of the alloy divided 
by the area of the strip or the equivalent cylindrical wire. 
Now substituting equation (6.2) into equation (6.2) yields: 



- E + Z alloy J 


( Z -Z , . ) 

mn alloy 



(6.3) 


or 


(6.4) 
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The conjugate gradient method can be employed next to solve 
for the currents as described in Chapter 4. 

Another objective of the study here is to relate the 
reflection coefficient evaluated by this method to the 
emmissivity and reflectivity measurements carried out by 
NASA, at the Langley Research Center. 

In conclusion, two techniques were developed here with 
an eye toward more efficient numerical computation for 
grating and mesh scattering. The first method, the Spectral 
Iteration Approach is applied to regions where the spacing 
between the wires is not less than two wavelengths. The 
second method, the Spectral Domain Conjugate Gradient 
Method, can be used for any spacing. Both techniques were 
solved in the Spectral Domain and both follow from one basic 
formulation. A study of the electromagnetic properties such 
as reflection coefficients, induced currents and aperture 
fields were presented and compared with data calculated by 
other methods to support the validity of the algorithm. 

A number of suggestions for possible extensions of the 
current algorithm to solve the problems of skew-symmetric 
structures, double screens, wires made of alloys with 
different resistivities and a mesh above ground were 
mentioned. The code used in the Fortran program and a 
listings of all main and utility subroutines appear in the 
Appendices . 
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8 . APPENDICES 

8.1 DERIVATION OF THE EQUATION FOR H AS A FUNCTION OF THE 
ELECTRIC VECTOR POTENTIAL F 

This appendix is to derive equation (2.9) from equation 
(2.7) in Chapter 2. 

Start with the equation: 

1 

E = - — V x F (8.1.1) 

Substituting equation (8.1.1) into the following Maxwell's 
equation: 

V X H = j toe E . (8.1.2) 

yields : 

V x H = -j cdVx F * (8.1.3) 

Equation (8.1.3) can be written as: 

y x (H + j <o?) = 0 (8.1.4) 

Now using the vector identity 

=0 (8.1.5) 

equation (8.1.4) becomes: 

H = -j co F -V® (8.1.6) 

v m 

where is the magnetic scalar potential. Taking the 

curl of equation (8.1.1) leads to: 

VxE V*V * F (8.1.7) 

which can be written as: 

VxE = -— tV V* ^ -V 2 (8.1.8) 
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by making use of the vector identity 

y x \7 x ! = y < y. ? ) - y 2 f (8.1.9) 

To completely specify the vector F, its divergence and 
its curl must be defined. In equation (8.1.1) the curl of F 
was defined. Now, one is at liberty to define the 
divergence of "f, which is independent of its curl. The 
choice of ^.F^ is made to simplify equation (8.1.8) which is 
achieved by letting: 

SJ • F = - j cjeu ® m (8.1.10) 

which gives: 


© 


m 


V-F 


( 8 . 1 . 11 ) 


3 uue 

Substituting equation (8.1.11) into equation (8.1.6) leads 
to: 


H - -j«P - V ( V • F) 

j <oue 
or 

H = -jcoF + \ 7 V- F 

j coue 

which is the same as equation (2.9) in Chapter 2. 
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8.2 DERIVATION OF THE EQUATION FOR THE SCATTERED MAGNETIC 
FIELD H S 

The purpose of this appendix is to derive equation 
(2.12-b) from (2.12-a) in Chapter 2. 

Start from equation (2.12-a), i.e. 


H S+ =- 


}cou 



k 2 - n 2 

*o a mn 


a 6 
mn M mn 


" a mn ^ mn 

I - __ 

G (a.__ ) M(a mri , S mn ) 

mn mn mn mn 

.2 2 

K o " mn , , . , L _ . , 

_J - ex P[3 ( 


( 8 . 2 . 1 ) 


Figure (8.2.1) below shows how the equivalence theorem could 
be utilized to transform the free standing inductive surface 
into a perfect electrical conductor [16] . 
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4 



ape.|cond.| 
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for z^O 




perfect electrical 
conductor 



image of 
E 1 due to 
electr. cond. 



-fee- 



-444- X 


Fig . 8.2.1 


Equivalent problem for an inductive FSS 
structure . 
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For the region z>0 the total H field (H. . ) at z=0 can be 

tot 


expressed as: 


H 


rs+ 


tot = H u 'y> + H t 


me 


( 8 . 2 . 2 ) 


:a .. 


Now the magnetic current related to the aperture field E is 
given by: 


■* + "Sa a 
M = E x n 


(8.2.3) 


, A A A 

where n is the normal to the aperture. For z > 0 n = z 

, - A A 

and for z < 0 n = -z. So 


A A . a a a A. A, 

M x x + My y = [E x x + E y y] x (z) 


a a a A 
= -E x y + E y x 


Similarly for z < 0 
M = E a x (-z) 


(8.2.4) 


(8.2.5) 


and the total H field is given by: 
-* s _ 

H tot - H 


( 8 . 2 . 6 ) 


At 2 " 0 H tot = H tot 


or 


A 3 ' = H S+ + H. lnC 


(8.2.7) 


Moreover -H S ~ = H S+ 


( 8 . 2 . 8 ) 


So equation (8.2.7) becomes: 
-2 H S+ = H t inC 


(8.2.9) 
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From the previous figures we had 2 M (Due to image theory) 
So multiplying equation (8.2.1) by a factor of two yields: 


- 2 i s+ =- L'Y' 

joou-Z I 


k 2 - -a 3 

o mn mn mn 


-ct mn®mn **o ® mn 


G M + exp[j(a mn x+6 mn y)l 


or 


H, 


me 


JWU 


. 2 2 
^o ^n 


■a 3 
mn mn 


-a 3 
mn mn 


k 2 — a 2 
o ^ mn 


M. 


M 


exp [ j (a x+3 y ) ] 
mn mm 


( 8 . 2 . 11 ) 


But the equation (8.2.4) M = E and M = -E so equation 

x y y x 

(8.2.11) leads to: 


v nc =— y 


k - a - a 3 
o mn mn k mn 


2 2 

h a mn ^mn k o ^ mn 


-E. 


exp[j<a mn x +e nl nl ' )1 


( 8 . 2 . 12 ) 


or 


H. 


me 


=--7 

3 a)u/__ 


-a 3 k 2 - a 2 
Tnn mn o mn 


2 2 

k “S3 - a 3 

o mn mn K mn 


-E. 


ex P { j |a .n rt M l' 11 


(8.2.13) 


which is the same as equation (2.12) in Chapter 2. 
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8.3 MINIMIZATION IN THE DOMAIN FOR THE CONJUGATE GRADIENT 
METHOD 


The algorithm that minimizes the error functional ERRF = 


* * —] 
r (AA ) r 


1/2 4 


is the following: 


V 0) = z J (0) + E 1 

ran 


for n=0 


p(°) = Z* r (0) 

* mn 

The equations for the n th iteration are: 


a = 
n 



£<n) 

to 


-<n) 

2 


- (n+1) = ?< n) + a £ (n) 


n 


"£( n+1 ) 


^n " 


" U) ' <*n Z mn 



"J?( n+1 ) 

2 


"r ( n ) 

2 


p( n+1 ) 


= Z* ? (n+1) +. B n P (n) 


mn 


n = n+1 



8.4 CODE 



— i i ~i r~ 


Fig. 8.4.1 Square cell 


AA= Distance between centers of vertical strips in x 
direction (INPUT) 

BB= Distance between inner edges of vertical strips in x 
direction (INPUT) 

CC= Distance between centers of horizontal strips in y 
direction (INPUT) 

DD= Distance between inner edges of horizontal strips in 
direction (INPUT) 

F= Frequency (INPUT) ^ 


I0PT=1 

for 

rectangular meshes 

( INPUT) 

IOPT=0 

for 

parallel wire grids 

( INPUT) 

PSI 

is . 

angle Q 

( INPUT) 

ITM=0 

for 

TE polarization 

( INPUT) 

ITM=1 

for 

TM polarization 

( INPUT) 

NOI = 

Number 

of iterations 

( INPUT) 



IX = Sampling length for FFT (INPUT) 

ALAMB = Wavelength 

NX = number of samples incident in the aperture along x 
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Fig. 8.4.2. Sampling Arrangement 

NXl r NX2 = Edge points in x direction 
NY1,NY2 * Edge points in y direction 


V,U * Expressions of the 

Floquet modes a 

mn 

and 

^ mn 

-j " 

J k 2 -(V 2 +U 2 ) 

for k 2 

> U 2 

+ V 

G = 

/ k 2 - < V 2 +U 2 ) 

for k 2 

< u 2 

+ v 


EXI , EYI = x and y components of the incident electric field 
( INPUT) 

HXI , HYI = x and y components of the incident magnetic field 
(INPUT) 

FFT3D = 3 dimensional complex Fast Fourier transform 
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FOR SIT 

X,Y 2 dimensional arrays for the aperture field 

JCX,JCY = 2 dimensional arrays for the current densities 
CONX , CONY = 2 dimensional arrays to store the constant: 

C -F ” 1 Z" 1 F[-H t lnc + e ( H t xnc ) ] 

XIX,YIY,XIY,YIX arrays used to store the perturbed aperture 
fields. 

GX,GY,HX,HY partial derivatives used in the contraction 
operation 
FOR THE S.D.C.G. 

ZINT internal impedance 
DX , DY direction vectors' 

RX , RY residual vectors 

X,Y are either the aperture fields, or the induced 
currents (the unknown) 

TX,TY two dimensional arrays used to store different values 
AN .v 

BN B n 

PHASEX , PHASE Y arrays used to store phase information 
CREFX,CREFY x and y components of the reflection 

coefficient 

REFF, REFT, RETT, RETF reflection coefficients along theta and 

phi angles for TE TM polarizations 
Z1,Z2 one dimensional arrays used to store the amplitude of 
the unknown X and Y so that they can be used for 
plotting purposes for the cross-unit cell 
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f 



Fig. 8.4.3. Cross-Unit Cell 

MDOX left edge point of perpendicular strips 
MUPX right edge point of perpendicular strips 
MDOY lower edge point of horizontal strips 
upper edge point of horizontal strips 


MUPY 
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8.5 FLOW-CHART FOR THE SPECTRAL DOMAIN CONJUGATE GRADIENT 
METHOD 



0 


onnr(O) 



Evaluate ERRF (0) 



^ * 


Calculate 

See Chapter 3 & 4 


'l 

J (n+l) -J {n) + a ^p(n) or 

E { n+1 > _g{ n ) + a ^p t n ) 

Z V ' 

Find ERRF (n+1) 

(See Chapter 3 & 4) 




Define the 
( n+ 

vector r 

new residual 
1) 


Calculate t 

he factor ft 
^ n 


V 



YES 

V 


Calculate the amplitude 
and phase of current 
densities and electric 
aperture fields. 

— -g — 

Evaluate reflection 
Coefficients 
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8.6 LISTING OF THE S.I.T. METHOD 


C«****sl?.FORT***** 

COMPLEX COHE,CZEEO,CDET,CX1 M.COHS.CXtlM.CSEFX.CREPY 
COMPLEX CONX1(32,32)/1024*{0. 0,0.0)/ 

COMPLEX CONY1 (32 >32) /1024* (0.0, 0.0) / 

COMPLEX CONX2 (32,32) /1024* (0.0, 0.0) / 

COMPLEX CONY2(32,32) /1024* (0.0, 0.0) / 

COMPLEX COJIX (32,32) /1024* (0.0, 0.0) / 

COMPLEX CONY (32,32) /1 024* (0 .0, 0. 0) / 

COMPLEX XU(32,32)/1024*(0. 0,0.0) / 

COMPLEX YO (32, 32) / 1 024* (0.0, 0.0) / 

COMPLEX G (32, 32) /1024* ( 0. 0,0.0)/ 

COMPLEX JCX (32,32) ,JCY (32,32) 

COMPLEX Y (32,32) /1024* (0.0, 0.0) / 

COMPLEX X(32,32)/1024*(0. 0,0.0)/ 

COMPLEX YIX(32, 32) /1024*(0. 0,0.0) / 

COMPLEX XIX (32,32) /1024* (0.0, 0.0)/ 

COMPLEX YIY (32 ,32) /1 024* (0.0 ,0.0)/ 

COMPLEX XIY (32,32) /1024* (0.0, 0.0) / 

COMPLEX GX (32, 32) , GY (32, 32) ,HX (32,32) ,HY p2,32) 

COMPLEX J,HX1,HYI,CWK(32) , & 1 1, A1 2, A21 , A22, DENO 
HEAL K ,X2 , RWK (342) 

DIMENSION AHP(32) , HIND EX (32) , INK (342) , CEOS (32) 

BEAL O (32) /32*0. 0/ 

BEAL Y(32,32)/1024*0.0/ 

C *** AA-DISTANCE BETWEEN CENTERS OP 7ERTIC. STRIPS IN X-DIRECTION ** 
C *** BB=DISTANCE BETWEEN INNES EDGES OP 7SBTIC. ST SIPS IN X-DIBECT.** 
C *** CC*DISTASCE BETWEEN CENTERS OF HORIZ. STRIPS IN Y-DIRECTION •• 

C *** OO^DISTANCE BETWEEN INNER EDGES OP HORIZ. STRIPS IN X-DIBECT.** 

READ (1,22) AA,BB,CC,DD,F,ERR 
22 FORMAT (8E1 0.4) 

P*2. 998E*8 

C *** IOPT=0 POE A RECTANGULAR OR SQUARE MESH ***** 

C *** IO PT* 1 POE A PARR ALLEL GRID ***** 

IOPT* 1 

IP(IOPT.GT.O) CC»1.500E*15 
IP(IOPT.GT.O) DD*1.500E*15 
WRITE (3, 33) AA ,.B3, CC ,DD, ERB 
33 FORMAT (' 0* , ' A* *,715.8,* B= *,P15.8,» C* ',F15.8, 

3* D= * ,P15. 8, • ERR* *,P15.8) 

WRITE{3,44) P 

44 FORMAT ('O', * FREQ = *,E10.4) 

READ (1 ,22) PHI, THI, PSI 
WRITE (3,55) PHI, THI, PSI 

55 FORMA? (* 9* , ' PHI® ',P10. 1,* THETA® ',P10.1,» PSI® *,P10. 1) 

C *** READ THE NUMBER OP SAMPLING POINTS **** 

READ (1,66) IX 

C *** ITH* 1 FOR ?H POLARIZATION **•**•♦* 

READ (1,66) ITM 

C *** READ HUMBER OF ITERATIONS ******* 

RE AD ( 1 , 66) NO I 

56 PORM AT (13) 
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Bl-2 .1415^3 
PI2=?I/2. 

TPI=6. 283135 
CV=2 . 9 97956 E+9 
tjy=4.E-7*PI 
FTD=57. 29578 
EP=8.854E-12 
E?A=SQRT (OO/EP) 

J=CMPIX (0.0,1. 0) 

ITER =1 

CONE=CMPLX (1.0, 0.0) 

C2ER0=CMPLX (0. 0,0. 0) 

W=TPI*F 

A1AM8=CV/F 

AA= A A/ALAtlB 

BB=BB/ALAMB 

CC=CC/ALAMB 

DD=DD/ALAMB 

C *** DETERMINE SAMPLING POINTS THAT CORRESPOND TO THE 
C CONDUCTING REGIONS AND THE APERTORE ***•**•♦* 

NX=IFIX(BB/AA*?LOAT(lX) *2.)/4*2 
NY=IFIX (DD/CC*FLOAT(lX) *2.)/4*2 
NX1= (IX-NX) /2+ 1 
NX2=NXUNX-1 
NY1= (IX-NY)/2*1 
NY2-NY1*NY-1 

WRITE {3 , 1 00) NX, NX 1 , NX2 , NY , NY1 , NY2 
100 FO RH AT ( • 0 ' , * NX**, 13,3 X, * NX1 a *,l3,3X,*HX2*',I3#-3X, 

. • NY=* ,I3,3X,'NY1 =, ,I3,3X, , NY2='»I3) 

K=TPI/ALAMB 

K2»K**2 

STSPK= SIN (THI/RTD) *SIN (PHI/RTD) *K 
STCPK=SIN (THI/RTD) *COS (PHI/RTD) *K 
CPS=COS (PSI/RTD) /SIN (PSI/RTD) 

110 CONTINOE 

C ***** CALCOLATE FLOQOET MODES ****♦*•*♦ 

DO 200 M=1,IX 
IP (M.GT.IX/2+1) GOTO 125 
D (M) =TPI* (M-1) /AA-STCPK 
GOTO 127 

125 tJ(M) =TPI* (M-IX-1) /AA-STCPK 
127 CONTINOE 

DO 190 N=1 , IX 

IF (M.GT. IX/2*1.AND.N.GT.IX/2*1) GO TO 160 
IF (M.GT. IX/2*1) GO TO 150 
IF(N.GT. IX/2+1) GO TO 140 

7(H,N)=TPI* (N- 1) /CC-TPI* (H-1)/A A*CPS-STSPK 
GO TO 170 

140 V (M, N) =TPI* (N-IX-1) /CC-TPI* (M-1) /AA*CPS-STSPK 
GO TO 170 

1 50 V (M , II) =TPI* (N-1) /CC-TPI* (M-IX-1) /A A*CP5-STSPK 
GO TO 170 

160 7 (M , H) =TPI* (N-IX-1) /CC-TPI* (M-IX-1 ) /AA*CPS-STS PK 

170 IP (K2.GE.0 (M)**2*V (M, N) **2) G(N, N) =-J*SQRT (K2- (0 (M)**2*V (H,N) **2 
• ) ) 

IF (R2.LT.U (H) **2*7 (M,N) **2) G (M, N) =-SQRT (0 (M) **2*7 (M,N)**2-K2) 
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. tcnitn 

190 CONTINUE 
200 CONTI!! as 

IF (ITX.GT. 0) GO TO 210 

c **** INCIDENT fields FOR te folariz ATION **** 

EX 1= SIN (-PHI/RTD) 

EY I=COS (PHI/RTO) 

HXI=COS (PHI/ETD) *COS (THI/RTD) /ETA 
HYI=SIN (PHI/STD) *COS (THI/RTD) /ETA 
GO TO 261 

C ***** INCIDENT FIELDS FOR TN POLARIZATION ****** 

210 EXI=COS (PHI/RTD) *COS (THI/RTD) 

EYI=SIN (-PHI/RTD) *COS (THI/RT D) 

HYI=SIN (PHI/RTD* PI2) /ETA 
HXI=COS (PHI/P.TD-PI21 /ETA 
261 COHTINDE 

C *** GIVE A GUESS POR INITIAL APERTORE FIELDS X 6 I •** 

DO 310 H=NX1,NX2 
DO 300 N=NY1,NY2 
X(H,H)='BXI 
Y (H,N) *EYI 
XT! (n,N)=X(H,N) 

YU (M,N)=Y(H,N) 

300 CONTINUE 
310 CONTINUE 

c **** START THE COMPUTATION OF CONSTANTS CONX 6 CONY THAT DEPEND 
C ON A GIVEN INCIDENT FIELD ** 

DO 315 1=1, IX 
DO 315 L=1,IX 
CONX1 (I,L) **1. *HXI*H*UU/J 
CONY 1 (I, L) = ♦ 1 . *HYI*W*UU/J 
CONX2(I,L)=-1.*HXI*W*UU/J 
CO MY 2 (I,L) *- 1. *HYI*B*UU/J 
315 CONTINUE 

C *** PERFORM THE TRUNCATION OPERATION *** 

DO 325 I=NX 1 , NX2 
DO 325 L=NY1,NY2 
CONX1 (I,L)=CZERO 
325 CONY1 (I,L) =CZERO 

C *** TAKE THE FOURIER TRANSFORM **** 

CALL PFT3D(CONX1,IX,IX,IX,IX,1,69,IVK,B«K,CNK) 

CALL PFT3D(CONY1 ,IX, IX, IX, IX, 1 , 69 , IHK,BVK,CNK) 

CALL FFT3D(CONX2,IX,IX,IX,IX,1,69,IHK,BVK,CWK) 

CALL FFT3D (CON Y2 ,IX,IX,IX,IX,1,69,IWK,RHK,CNK) 

DO 340 H=1 , IX 
DO .-340 N=1 , IX 

CONX (N,N) =CONX1 (3,N) 4-CONX2 (M,N) 

340 CONY (H,H) =CONY 1 (H,N) ♦CONY2 (M,N) 

DO 350 *1=1 ,IX 
DO 350 N=1 , IX 

CDET=- (U (fl) *V (H,H) /G (M, N) J **2- (7 (H , N) **2/G (11, N) -G ( H, N) ) 

. *(G(H,N)-H (N)**2/G(’1,N)) 

CONS=CONX(M,N) 

CONX (N,N}= (-U(N) *V (.1 ,N) /G(M, N) *CONX (3,N) - (V (!!,N) **2/G(N,N) -G (N,N)) 
. *CON Y (S,N) ) /CDET 

CONY (M, N) = (- (G (M , N) -U (fl) **2/G (H, N) ) *CONS*U(fl) *V(f», N)/G (f!,N) * 
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•CONY (M, N) ) /CDET 
350 CONTINUE 

c **** hi; d of calculation of CONX 6 CONY ***** 
c 

c **«* now START ITERATIVE PROCESS !!!!!!!!!!!!!!! * ***** 

C 

C *** SET ALL PARTIAL DERIVATIVES EQUAL TO ZERO ***** 

000 DO 541 H=1,IX 
DO 541 N=1,IX 
XIX (!1,N) =CZEBO 
YIX (B,N) =CZERO 
XIY (B,H) =CZERO 
YIY (B,N) =CZSRO 
541 CONTINUE 

C **** P2RT0RB APERTURE FIELDS BY (0.01,0.01) ******* 

DO 543 B=NX1,NX2 
DO 543 N=>NY1,NY2 
XIX (fl,N)=X (B,N) *(0.010,0.010) 

YIX (H,N)=Y (B,N) 

XI Y (H,N)=X(B,H) 

YlY (H,N)=Y (B,N) *(0.010,0.010) 

543 CONTINUE 

C *** TARE THE FOURIER TRANSPOSE OP THE APERTURE FIELDS X 6 I **** 
CALL FFT3D(X,IX,IX,IX,IX,1,69,IWK,RWK,CWK) 

CALI PFT3D (Y, IX, IX, IX, IX ,1 ,69, IWR, RWF ,CWK) 
c **** rtHLTIPLY TRANFORHED FIELDS BY THE PLOQUET COEFFICIENTS 
C AND GREEN* S PUNCTION **** 

DO 560 B=1,IX 
DO 550 N=1,IX 
CXBN=X(H,N) 

X (B,N)= (U(B) * V ( H , N) /G (H,N) *X (B, W) ♦ (V (B,N) **2/G(B, N) -G (H,N) ) 

. *Y(«,N)) 

Y (B,N) = ( (G (B, N) -U (B) **2/G (B, N) ) *CXHS-U (H) *V (R,N) /G (B, N) 

. *Y(B,N)) 

550 COHTIHUE 

560 CONTINUE 

C ***** TAKE THE INVERSE FOURIER TRANFORN ****** 

CALL FFT3D (X, IX, IX, IX, IX , 1 , -69,1 WK ,RWK,CWK) 

CALL FFT3D (Y, IX, IX, IX, IX, 1,-69, I WK , RWK,CWK) 

WRITE (3,570) ITER 

570 FORHAT(3X,/ • ITERATION RUBBER ',12) 

C *** CALCULATE CURRENT DENSITIES ***** 

DO 600 fl= 1, IX 
DO 600 N=1,IX 

jcx (n, N) =(y (h ,n) *j/w/uu*h¥i|»(-2.) 

600 JCYfH, N)=fr (B,N) *J/W/UU*nXl)*(2J 
C *** PLOT CURENTS ON STRIPS ***** 

DO 620 1=1, IX 
ABP(I) =CABS (JCY (1,1) ) 

EINDEX (I) = (FLOAT (I-IX/2) -. 5) /IX*AA*1.045 
IF (ITER.GT. (NOI-1) ) WRITE(8,*) AHP (I) , RIND EX (I) 

620 CONTINUE 

C IF (ITER.GT. (NOI-1) ) CALL GEHPT (RINDEX, AMP, IX , 0) 

C **** TRUNCATION **** 

DO 740 fl=NX1,NX2 
DO 730 H=NY1 , NY2 
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X (M, N) =CZEP.O 

Y ( M , II ) =CZERO 
73 0 CONTINUE 
740 CONTINUE 

C **** NOB FIND THE P (Tfi.ONC (INVERSE P (G E) ) ) **** 

C 

• CALL FFT3D(X, IX, IX, IX, IX, 1,69, IWK, RWK, CWK) 

CALL FPT3D (Y , IX, IX, IX, IX , 1 , 6 9, IWK, R WK ,CWK) 

I?ER=ITER+1 
90 760 H= 1 , IX 
DO 750 N=1,IX 

CDET=- (U(H) *V(H,N) /3 (H,N)) **2- (V (H, N) **2/G (H , N) -G (H. N) ) 

. * (G (K, N) -0 (H) **2/0(11, N) ) 

CXMN=X(fl,N) 

X (H, N) = (“U (N) * V (H , N) /G(H,N) *X(H,N) -(V (H, N) **2/G(H,N) -G (H,N)) 

• *Y (H, NJ ) /CDET 

Y(H,N) = (-(G (H, N) -U ( H) **2/G(fl,N) ) *CXHH + U(H) *V (H,N)/G(H, N) * 

. Y(H,N))/CDET 
750 CONTINUE 
760 CONTINUE 

C **** ADD F (TSUN (INVERSE F (G E) ) ) TO CONX AND CONY **** 

DO 780 H=1,IX 
DO 770 N=1,IX 
X (N, N) =X (H, 5) ♦CONX (S,N) 

Y (H, N) =Y (H, N) ♦CONY (H,N) 

770 CONTINUE 

780 CONTINUE 

C *** CALCULATE THE REFLECTION COEFFICIENTS CREFX AND CREFY **♦ 
CRSFX=X (1 , 1 ) / (FLOAT (IX) * FLOAT (IX) ) -EXI 
CREFY® Y (1,1) / (FLOAT ( IX) *FLOAT(IX)) -EYI 
REFX=CABS (CREFX) 

RSFY=CABS (CREFY) 

WRITE(3, 800) REFY 
800 FORHAT (10X,2P10.3) 

C *** TAKE THE INVERSE FOURIER TBANSFORH OF THE RESULT TO OBTAIN 
C A NEW VALTTE FOR THE APERTURE FIELDS **** 

CALL FFT3D ( X,IX , IX , IX , IX , 1 , -69 , I WK , RWK , CVK) 

CALL FFT3D (Y, IX, IX, IX, IX ,1,-69, IWK, RWK, CWK) 

C *** PLOT APERTURE FIELD *** 

DO 830 1=1, IX 

AHP (I) =CABS (Y(I, 16) ) 

C CROS (I) =CABS (Y (8,1) ) 

RINDEX (I) = (FLOAT (I-IX/2) 5) /IX* AA*1. 045 

IF(ITER.EQ.NOI) WRITE(8,*) AHP (I) , RINDEX (I) 

83 0 CONTINUE 

C IF (ITER. GE. NOI) CALL GENPT (RINDEX, AHP, IX, 0) 

C CALL GENPT (RINDEX, CROS, IX, 0) 

C *** REPEAT SAHE PROCESS FOR PERTURBED FIELDS XIX ,XIY,YIX,S YIY *** 
CALL FFT3D (XIX, IX, IX, IX, IX, 1 ,69 ,IWK,RWK, CWK) 

CALL FFT3D (YIX, IX, IX, IX, IX, 1,69, IWK, RWK, CWK) 

CALL FFT3D ( XI Y, IX, IX, IX, IX, 1,69, IWK, RWK, CWK) 

CALL FFT3D (YIY, IX, IX, IX, IX, 1 , 69, IWK,Rf K,CWK) 

DO 850 H= 1 , IX 
DO 840 N= 1 , IX 
CXHN=XIX(H,N) 

XIX (H,N)= (U(.H) *V (M,N) /G <H, N) *XIX <M,N) ♦ (V (M, N) **2/G (H, N) -3 (H, N) ) 
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. *UX(n,N)) 

Y IX (B,M) = ( (G(M, N) -U (M) **2/G ) *CXBtWI (N) *V (M,N) /G(B,N) 

. *YIX(B,H)) 

CXIMN=XIY (B,N) 

XI Y (B,N) = (0 (.1) * V (B, N) /G (,*1, N) *XIY (M, N) ♦ (V (.1, N) **2/G(M, N) -G (B,N) ) 
. *YIY(M,N)} 

YIY (M,tJ) = { (G (B, N) -U (1) **2/G (B , N) ) *CXIMN-U (B) * V (B, N) /G (B, N) 

. *YIY(B,N)) 

940 CONTINUE 

950 CONTINUE 

CALI FFT3D{XIX # IX,IX,IX,IX,1,-69,IWK,RWK,CVK) 

CALL FFT3D(YIX,IX,IX,IX,IX,1,-69,IHK,RWK,CWK) 

CALL FFT3D{XIY,IX,IX,IX # IX,1,-69,IWK,RHK,CWK) 

CALL FPT3D (YIY, IX, IX, IX, IX, 1 «~69 ,1 WK, RWK.CWK) 

,C 

C*** N OS PERFORM THE TRUNCATION OPERATION FOR PERTURBED FIELDS *** 

C 

DO 870 H=NX1,NX2 
DO 860 N-NY1 , NY2 
XIX (S,N)=C2ERO 
YIX (H,N) =C2ER0 
XIY {H/ N) =C2ERO 
YIY (B,N)=C2ERO 
860 CONTINUE 
870 CONTINUE 
C 

c **** HOW FIND THE P (TRONC (INVERSE F (G EP) ) ) *** 

C 

CALL PPT3D(XIX f IX,IX # IX,IX,1,69,IWK,RWK,CWK) 

CALL FFT3D (YIX,IX,IX,IX,IX,1 ,69,IWK,RHK,CHX) 

- CALL PFT3D (XIY , IX, IX, IX, IX , 1 ,69,IWK,RWK,CWK) 

CALL PFT3D(IIY,IX,IX,IX,IX,1,69,IWK,BWK,CWK) 

DO 910 B 3 1 fIX 
DO 900 N* 1 , IX 

CDET=- (U(B)*V(fl,N) /G ( B , N) ) **2- (V (S,N)**2/G (B,N)-G(H,H) ) 

. * (G (H, N) -0 (H) **2/G (S,N) ) 

CXBN»XIX(f! # N) 

XIX («, N) « (-0 (B) *V (B, N) /G (B, N) *XIX (S, N) - (V (B, N) **2/G (B, N) -G (H,N) ) 
* *YIX (H,N))/CDET 

YIX (B, N) 3 (- (G (B,R) -0 (N) **2/G (B, N) ) *CXNN*0(H) *V (8,N)/G(H, N) * 

. YIX (B, N) ) /CDET 
CXIBH=XIY (B,N) 

XIY (fl,N) = (-U(B) *V(M,N) /G (B, N) *XIY (H,N) - (V <B, N) **2/G (B, N) -G (B, N) ) 
.*YIY(H,N))/CDET 

YIY (B, N) s {- (G (B, R) -U (H) **2/G (B, N) ) *CXIHN*U (BJ *V (B, N) /G (H, N) * 

. YIY (H,N)) /CDET 
900 CONTINUE 
910 CONTINUE 

DO 930 B*1,IX 
DO 920 N*1,IX 

XIX (8# N) =XIX (B,N) +CONX (B,N) 

YIX (B, N)=YIX (M,H) «-CONY (B, N) 

XIY (H,N) =XIY(fl,N) +COMX (B,N) 

YIY (B, N) =YI Y (B ,H) ♦CONY (B, H) 

920 CONTINUE 
930 CONTINUE 



ORIGINAL PAGE IS 
OF PbOR QUALITY 


CALL FFT3D (XIX, IX, IX, IX, IX , 1 , -6 9,1 WK, RKK,C WK) 

CALL FFT3D(YIX,IX,IX,IX,IX,1,-69,IWK,BWK,CWK) 

CALL FFT3D(XIY,IX,IX,IX,IX,1 ,-69,IWK,RWK,CWK) 

CALL FFT3D(YIY,IX,IX,IX,IX,1,-69,IWK,RWK,CWK) 

C *** EVALUATE PARTIAL DERIVATIVES GX,GY,HX,SHY **** 

DO 950 5=1, IX 
DO 940 N= 1 ,IX 

GX (M, fJ) = (XIX (H, N) -X (M, N) ) / (0. 0 10 ,0. 0 10) 

HX (R,N)=(YIX(R,N)-Y(R,N) ) / (0. 01 0 , 0 . 0 10) 

GY (fl, N) = (XIY(R,N)-X(8, N) ) / (0. 0 1 0, 0. 01 0) 

94 0 HY (H, N) - (YIY (R , N) - Y (8, N) ) / (0. 0 1 0 ,0. 0 10) 

950 CONTINUE 

C *+* IMPROVE PREVIOUS ITERATE FOR APERTURE PIELDS BY USING 
C A CONTRACTION FACTOR **** 

DO 960 R*NX1 # NX2 
DO 960 N=NY1 ,NY2 

DEKO=GX (R,N) *HY (H, N) -HX (8, N) *GY (R, N) -GX (H,N) -HY(H,N) ♦!. 

A 11= (1.-HY (R, N) ) /DENO 
A12=-GY(R,N) /DENO 
A21=-HX(H,N) /DENO 
A22= (1.-GX (R, N) ) /DENO 
CXHN=X(R,N) 

X (R, N) =A1 1*X (R,N) *A 1-2*Y ( R* N) ♦ (1.-A11)*XO (H,N)-A12*Y0 (M,N) 

Y (R, N) =A2 1 *CXRN+ A22*Y (R, N) -A21*X0 (H, N) ♦ ( 1. -A22) * YO (H,N) 

960 CONTINUE 

IF (ITER. GT. NOI) GO TO 1000 
C . 

C ROW TRUNCATE THE IRPROVED APERTURE FIELD E ***» 

C 

DO 980 5=1, IX 
DO 970 N= 1 , IX 

IF(R.GE,NX1.AND.R.LE.NX2.AND.N.GE.NY1.AND.N.LB.NY2) GO TO 965 
X(fl,N) =(0.00,0.00) 

Y (R,N)= (0.00,0.00) 

965 XU (H ,N) = X (H , N) 

YU(R,N)=Y(H,H) 

970 CONTINUE 
980 CONTINUE 
GO TO 400 
1000 CONTINUE 
STOP 
END 
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8.7 LISTING OF THE S.D.C.G. METHOD FOR THIN STRIPS WITH 
THE SQUARE-SHAPED UNIT CELL 


C***»* THIS IS TIIK COWJG. GRAD. METHOD POR CIIP.R ENTS OH THIN WIRES** 
C**** MINIMIZATION IN THE RANGE (VAN DER BERG) ****** 

COMPLEX CO!lE,CZEPO,CXHN, FI 0 
COMPLEX CREFX, CREPT, CREF,CRET, ZINT 
COMPLEX G (32, 32)/1024*{0. 0,0.0)/ 

COMPLEX Y(32,32)/1024* (0. 0,0.0)/ 

COMPLEX X(32,32) /1024* (0.0, 0.0) / 

COMPLEX YU { 32, 32) /I 024 * ( 0. 0, 0. 0) / 

COMPLEX XU (32, 3 2) /10 24 *(0.0, 0.0) / 

COMPLEX RX(32,32)/1024*(0. 0,0.0) / 

COMPLEX RY(32, 32)/1024*(0. 0,0.0) / 

COMPLEX J, HXI, HYI,CW fC(32) 

COMPLEX DY (32,32) /1 024* (0.0, 0.0) / 

COMPLEX DX(32,32)/1024*(0. 0,0.0) / 

COMPLEX TX (32, 32)/1024*(0. 0,0.0) / 

COMPLEX TY (32, 32) /1024* (0.0, 0.0) / 

REAL K,K2, RWK ( 342) 

DIMENSION IWK(342) ,RR(3S0) ,CH{350) , PHASEX (32 ,3 2) ,PHASEY (32 ,32) , 
.21 (1024) ,X1 (1024) ,Y1 (1024) ,Z2 (1024) , AMP(36) , RINDEX (36) 

HEAL O (32) /32* 0.0/ 

REAL V (32,32)/1024*0.0/ 

C *** AA=DISTANCE BETWEEN CENTERS OF VBRTIC. STRIPS IN X-DIRBCTION *• 
C *** B3=DISTANCE BETWEEN INNER EDGES OF YERTIC. STRIPS IN X-DIRECT.** 

C *** CC=DISTANCE BETWEEN CENTERS OF EORIZ. STRIPS IN Y-DIHECTION ** 

RE AD ( 1 , 10) A A, BB,CC,DP,F, ERP 
10 ' FORMAT (8E10.4) 

... F=2 . 94 8E*8 

C *** IOPT=0 FOR A SQUARE OB A RECTI NGOLAR MESH **** 

C *** topt=1 FOR A PARALLEL GRID **** 

IOPT=1 

IF (IOPT.GT.0) CC= 1. 500E +1 5 
IF (IOPT.GT. 0) DD= 1 . 500E+1 5 
WRITE (3, 20) AA,BB,CC,DD,ERE 
20 FORMAT (' 0* , * A* «,F15.8,» B* ',715.8,* C= ',F15.8, 

.' D= ',F15.8,« ESP= ' , F15. 8) 

WRITE (3, 30) F 

30 FORMAT ( ' 0 * , ' FREQ = ',E10.4) 

READ (1 , 1 0) PHI ,THI , PSI 
WRITE (3, 40) PHI, THI, PSI 

40 FORM AT ( ' 0 ' , ' PHI* ',F10.1,» TIIETA* • ,P10. 1 ,» PSI= *,F10.1) 

C *** READ THE NUMBER OF SAMPLING POINTS ****** 

READ ( 1 , 50) IX 

C *** ITM=0 FOR TE POLARIZATION **** 

C *** ITM= 1 FOR TM POLARIZATION **** 

READ (1,50) ITM 
WRITE (3,45) ITM 

C *** READ THE NUMBER OF ITERATIONS *** 

READ (1,50) NOI 

45 FORMAT (IX, 'THE YALUE FOR ITM 15= ',.13) 

50 FORM AT (1 3) 

PI=3. 141593 
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?:.•? = r> r/2 . 

TPT=f>. 28 31 95 
CV=2.997956E*8 
rjrj=a.E-7*Pl 
RTD=57. 29578 
E?=8.8S4E- 12 
ETA=SQRT(nO/EP) 

J=CMPLX (0.0, 1.0) 

ITEF.= 1 

CONE=CHPLX (1. 0,0.0) 

CZEPO=CMPLX(0. 0,0.0) 

SIG? A=5. E20 
W=TPI*F 

ZINT= (1. 0, 1. 0) *SQP.T(W*UTJ/2. /SIGMA) /(1.0) 

ALAWB=CV/P 
AA= AA/ALAM B 
BB=BB/ALA«B 
CC=CC/ALAHB 
DD=DD/ALAHB 

NX=IPIX (BB/AA* FLOAT (IX) *2. ) /4*2 

NY=IFIX(DD/CC*PLOAT(IX) *2. )/«*2 

NX1= (I X-.NX) /2* 1 

NX2=NXHNX-1 

NY 1= (X X-NY) /2+ 1 

NY2=NYUNY-1 

WRITE (3, 60) NX,NX1,NX2,NY,NY1,NY2 
60 FORMAT ('O' , • NX*',I3,3X,'NX1=*,I3,3X,»NX2»* ,I3,3X, 

. • NT=* ,13, 3X , * NI1 = * , 13 ,3X, • NY2*' ,13) 

K=TPI/ALAHB 

K2=K**2 

STSPK= SIN (TIII/RTD) *SIN (PHI/RTD) *K 
STCPK=SIN (THI/RTD) *COS (PHI/RTD) *K 
C?S=COS (PSI/ETD) /SIN (PSI/RTD) 

70 CONTINUE 

C *** DEFINE THE FLOQOET COEFFICIENTS **** 

DO 100 H=1,IX 
IF (M.GT. IX/2+1 ) GO TO 75 
0(H) =T?I*(H-1) /AA-STCPK 
GO TO 80 

75 ri(M) =TPI* (K-IX-1) /AA-STCPK 

80 CONTINOE 

DO 90 N= 1 , IX 

IF (H.GT. IX/2*1 .AND.N.GT. IX/2*1) GO TO 84 
IF (K.GT. IX/2+1) GO TO 83 
IF (N.GT. IX/2*1) GO TO 81 

7 (M,N) =TPI* (N-1) /CC-TPI* (H-1) A A*CPS-STSPK 
GO TO 85 

81 V (W,N) =TPI* (N-IX-1) /CC-TPI* (H-1) /A A*CPS-STSPK 
GO TO 85 

83 V (H,N) =TPI* (H-1) /CC-TPI* (H-IX-1) /A A*CPS-STSPK 
GO TO 85 

84 V (N,N) =TPI* (N-IX-1) /CC-TPI* (H-IX-1) /AA*CPS-STSPK 

85 IF (K2 . GE. 0 (H) **2*7 (H,H) **2) G(H,N) =-J*SQRT(K2- (O (H) **2 *7 (H, N) **2 

.)) 

IF (E 2. LT.O («) **2*7 (M,N) **2) G (H, N) *-SQPT(P (M ) * *2+V (fj , N) **2-K2) 

. ‘CONE 
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90 COVTINUE 
100 CONTINUE 

IF (ITM.GT.O) GO TO 110 
C 

C *** INPUT FOR TE POLARIZATION **** 

FXI=SIM (-PHI/ETD) 

PYT= COS (PHI/RTD) 

HXI=COS (PHI/RTD) *COS (THI/RTD) /ET A 
HYI=SIN (PHI/RTD) *COS (TIII/RTD) /ETA 
EF=1.0 
GOTO 120 

C **** INPUT FOR TN POLARIZATION **** 

C 

110 EX I=CO S (PHI/RTD) *COS (THI/RTD) 

FYI=SIN (-PHI/RTD) *COS (THI/RTD) 

HYI=SIN( PHI/RTD- PI2) /ETA 
HXI=COS (PHI/RTD- PI 2) /ETA 
ET=1.0*COS (THI/PTD) 

120 CONTINUE 

c **** set YOUR INITIAL GUESS POR X AND Y **** 

C 

123 DO 130 H=1,IX 
DO 125 N=1,IX 
X (N, N) =CZEP.O 
Y {(1, N) =CZERO 
XU(H,N) = X{H,H) 

YU (H,N)=*Y(H,N) 

125 CONTINUE' 

130 CONTINUE 

c **** WORK ON INITIAL GUESS*** 

C 

CALL FPT3D (XU, IX, IX, IX, IX, 1,69,IWK,RWK,CWK) 

CALL PFT3D (YU, IX, IX, IX, IX, 1 , 69,IWK,RWK,CKK) 

DO 160 f!= 1 , IX 
DO 150 N=1,IX 
CXHN=XU(H,N) 

XU(N,N) = ( (G(N,N) -V (M, N) **2/G(N,N) )*XU(H,N)- (U (11) *V (H, N) /G (H,N)> 

. *YU(a,N))/(J*H*EP)/2 

YU (N, N) = (-U (N) *V (H,N) /G (H,N) *CXHN* (G (H, N) -U («) **2/G (H, N) ) 

. *YU(H,N))/(J*H*EP)/2. 

150 CONTINUE 
160 CONTINUE 

CALL FPT3D (XU, IX ,IX,IX ,IX, 1 , -69, INK, RHK, CWK) 

CALL FFT3D (YU, IX, IX, IX, IX, 1,-69, INK, RWK,CWK) 

C WRITE ( 3, 1 TO) ITER 

170 FORRAT(3X,/ • ITERATION NUNBER ',12) 

DO 200 H=1,IX 
DO 190 N=1 , IX 

C ** CONFUTE THE ERROR FOR YOUR INITIAL GUESS ** * 

RX (K,N)=EXI*XU (H,N) 

RY (fl, N) = EYT*YU (H,N) 

IP (M.GE.NX1 . AND.H. LE . NX2 . A »D . N.GE . NY 1 . AND. N.I.E.NY2) PX (S,N) =CZERO 
IF (N . GE. NX 1. AN D. N. LE* NX2 . A NO. N. GE. NY1 • AND. N« LE. NY2 ) RY (H,N) =CZ ERO 
ERROR= ERROR* RX (H, N) *CONJG(RX (N , N) ) *RY (N, N) *CONJG (RY (t!,N) ) 

F5 = F5*RX (M , N) *CONJG (RX (H , H) ) *R Y ( PI , N) *CONJG (RY (M, N) ) 

DX <H,N) =RX(fl,N' 
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DY (*,N) =PY (M,N) 

190 CONTINUE 
200 CONTINUE 
C 

•C «■*** FIND THE FOUP.IES TRANSFORM OF RESIDUAL *** 

C 

CALL FFT3D (DX, IX, IE, IX, IX, 1 , 69,IWK,RWK,CWK) 

CALL FFT3D ( DY, IX , IX, IX , IX, 1 , 6 9, I WK ,R WK,C WK) 

C ** flUTILPY BY THE CONJUGATE TRANSPOSE OF THE MATRIX Z 
C TO FIND THE VECTORS DX AND DY *** 

DO 220 M=1,IX 
DO 210 N=1 , IX 
CXNN=DX (M,N) 

D X ( M , N) ={ CO NJ G ( G ( M , N) -V (N, N) **2/G (H,N) ) *DX(M, N) - 
. CON.1G (V{M,N)*a(M) /G (M, N) ) *01 (1, N) ) /CONJG ( J* B*EP)/2. 

DY (M,N) = (CONJG(-V {fl # K) *U(M)/G (M,N) ) *CXHN* 

. CONJG (G (N,N) -U(M) **2/G (H,N) ) *DY (M,N)) /CONJG (J*W*EP)/2. 

210 CONTINUE 
220 CONTINUE 

C **** NOW FIND THE INV. FOURIER TEAS. OF THE DIREC. FUNCTIONS ** 

C 

CALL FFT3D(DX,IX,IX,IX,IX,1,-69,IWK,RWK,CWK) 

CALL FFT3D (DY,IX,IX, IX, IX ,1,“69,IWK,RWK,CWK) 

DO 360 H»1,IX 
DO 350 N=1, IX 

DX (N,N)=DX (H,N)-RX (M,N) *CONJG(ZINT) 

DY (M,N) = DY{H,N)-RY(M,N) *CONJG(ZINT) 

IF (M . GE. NX1 . AND. (1. LE. NX2. AND. N. GE. NY 1 . AND. N. LE.NY2 ) DX (M ,N) *CZERO 
IF (M.GE. NX1 . AND. N. LE. NX2 . A ND. N. GE. NY 1 . AND. N. LE.NY2) DY(f!,N) =CZERO 
TY (N,N) =DY (M,N) 

TX (M,N) =DX(M,N) 

F3=F 3+CONJG (DX (M,N)) *DX(M,N) *C0NJG (DY (M, N) )*DY(M,N) 

350 CONTINUE 
360 CONTINUE 

c **** now START THE ITERATION PROCESS (MINIMIZATION) *** 

C 

C **** MULTIPLY THE DIRECTION VECTOR BY THE MATRIX Z **** 

C **** STORE YOUR RESUTLS IN VECTORS TX AND TY **** 

365 CALL FFT3D(TX,IX,IX,IX,IX,1,69,IVK,RWK,CWK) 

CALL FFT3D (TY, IX, IX, IX, IX, 1, 69,IWK,RWK,CWK) 

DO 400 H=1,IX 
DO 370 N*1,IX 
CXRN=TX (M , N) 

TX(M,N) = ((G(N,N)-V(M,N) ** 2/G (M, N) 1 *TX (K,N)-<U (M)*V(H,N) /G (H, N) ) 

. •TY(M,N))/(J*W*EP)/2. 

TY (M, N) = (-U (M) *V (M,N) /G (M, N) *CXKN ♦ (G (M, H) -0 (S) **2/G (M,N) ) 

. *TY (M,N) )/(J*W*EP)/l. 

370 CONTINUE 

400 CONTINUE 

CALL PFT3D(TX,IX,IX,IX,IX, 1, -69 ,IVFC, RWK,CWK) 

CALL FFT3D (TY, IX ,T X, IT , IX, 1 , -6 9 , TU K, EWK, CWK) 

F1=0. 0 

DO 410 M=1 , IX 
DO 410 N=1 , IX 

TX (fl,N) =TX (M,N) -DX (H, N) *ZINT 
TY (d,N) =TY (M#N) -DY (11, N) e ZINT 
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IP (fl.GB.MX1 . AN 0. fl. LE. NX2 . A ND . N. GE. NY 1 . AND. N.L2.NY2) TX (M,N)=CZER0 
IF (H.GE. NX1. AND. M. LE.NX2.AND.II.GE.MY1. AND. N. LS.MY2) TY (N,N) =CZER0 
410 FI =F1 ♦CONJG (TX (M,N) ) *TX (H,N) ♦CONJG (TY (M, N) )*TY (M,N) 

ITEEaITER+1 

C **** COMPUTE CONSTANT AN •*** 

an=f.vfi 

CH (ITER) =SQET( ERROR) /SQRT (FS) 

ESR=CH (ITES) *1 00 
C •** CALCULATE ERSOB **** 

ERHOR= ERROR- (P3**2/F1) 

C *•* GET A NEW ESTIMATE FOR YOUR UNKNOWNS X C Y **** 

DO 560 H=1 # IX 
DO 550 N=1,IX 
X (M,N) =X (M,N) ♦ AN*DX(M,N) 

Y(M,N)=Y(M,N)*AN*DY(M, N) 

550 CONTINUE 
560 CONTINUE 

C **** A NEW ESTIMATE FOR THE RESIDUAL VECTORS RX S BY *•** 

DO 580 H=1, IX 
DO 570 N=*1,IX 
RX (H,N)aRX (H,N)-AN*TX(M,N) 

RY(H,N)aRY (H,N)-AN*TY(M,N) 

TX (M,N) =RX (M,N) 

TY (H ,N) aRY (H,N) 

570 CONTINUE 
580 CONTINUE 

HR (ITER) aptOAT (ITER) 

WRITE (8 ,*) CH (ITEP) ,RR (ITER) 

C **** MULTIPLY THE RESIDUAL VECTORS BY THE CONJG. TRANS. OF MATRIX Z •* 
CALL FPT3U<TX,IX,IX,IX,IX,1,69,IWK,RWK,CWK) 

CALL PPT3D (TY, IX, IX, IX, IX, 1, 69,IWK,RWK,CWK) 

DO 600 H«1 , IX 
DO 590 Nal.IX 
CXMNaTX (M,N) 

TX (M, N) a (CONJG (G (M, N) -V (M,N) **2/G (M,N) ) *TX (H, II) - 
. CONJG (V{M,N)*U(M) /G (M,N) ) *TY (M, N) ) /CONJG (J*W*EP)/3. 

TY (M# N) a (CONJG (-V (N, N) *U (H) /G (M,N) ) *CXMN* 

. CONJG (G (M, N) -U (H) **2/G (M, N) ) *TY (R,N) ) /CONJG (J*W*BP)/2. 

590 CONTINUE 
600 CONTINUE 

CALL FFT3D (TX,XX,IX,IX,IX, 1, -69,IWK,RWK, CWK) 

CALL FFT3D (TY, IX, IX, IX, IX, 1,-69 ,IWK,RVK, CWK) 

C ***** STORE THE OLD VALUE FOR 73 IN 72 TO CALULATE DN LATER **** 
P2*P3 
F3»0.0 

DO 644 Nal.IX 
DO 644 ?J=1 f IX 

TX (N,N) ®TX (M,N) -RX (M, N) *CONJG(ZINT) 

TY (M,N)aTY (M,N) -RY (M,N) *CONJG(ZINT) 

IF (M.GE. NX1. AND. M. LE. NX2 . A ND. N. GE. NY 1 . AND. N. LE. NY2) TX (M,N)aCZERO 
IF (M.GE. MX 1 . AND. M. LF. NX2. AND. N.GE. NY1 . AND. N. LE.NY2) TY (N,N)*CZEBO 
F3=F3*CONJG (TX (M , V) ) *TX (M,N) ♦CONJG (TY (M, N) )*TY (M,N) 

644 CONTINUE 
C NOW CALULATE 3N 

9N=E3/P2 


•**«*#* 
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Cl IF(ITEn.EQ.20.OR.ITER.EQ.40.OR.ITER.EQ.60) BN=0.0 

C *»** OBT AIN A N EV ESTIMATE FOR THP DIRECTION" VECTORS DX f, DY * 

00 664 N=1,IX 
no 654 N=1,IX 
PX(M f N)=?X(M,N) ♦BN#DX (M, N) 

OY (H,N)=TY,(M,N) ♦ BN*DY(!1, N) 

TX (M, N) = DX (3 ,M) 

TY (M,N) =DY (M,N) 

654 CONTI N BE 
664 CONTINUE 
C 

c **** CONTINUE THE ITERATIVE PROCESS *<■** 

IF {ITER. GT. NOT) CALL GENPT (RR,CH,NOI,0) 

IF (ITER. GT.NOI)GO TO 700 
IF (ERR.LT. 1) GO TO 700 
GO TO 365 

C ** STORE X SY INTO 1-DIHEN. ARRAYS TO BE USED FOR ANY PLOTTING 
C PURPOSES * 

700 DO 720 M= 1 ,IX 
DO 720 N=1,IX 
I=(H-1)*IX*N 

21 (I) =CABS(X(M,N)) 

22 (I) -CABS (Y (H ,N) ) 

XI II) = (FLOAT (H-IX/2) 5) /IX*AA*1.05 
Y1 (I) = (FLOAT(N-IX/2)-.5) /tX*BB*1. 05 
WRITE (7,*) XI (I) «Y1 (I) »Z1 (I) ,Z2 (I) 

720 CONTINUE 

GO TO 730 

C *** FIND THE PHASE FOR THE CURRENTS X & Y ***** 

723 DO .725 M= 1 , IX 
DO 725 N=1,IX 
PHASEX (M,N) =0.0 
PHASE! (N,N) =0. 0 

IF (H. GE. NX1 . AND. H. LE.NX2 . A ND. N. GE. NY1 . AND. N. LE. NY2) GO TO 725 
REX=REAL (X (H # N) ) 

AIHX= AIM AG (X (M, N) ) 

IF (REX.GE. 0. 0. AND. AIHX.GE. 0. ) PX=ATAN (AIHX/REX) *RTD 
IF <REX.LT. 0.0. AND. AIHX.GE. 0. ) PX= 180.-ATAN (AIMX/REX) •RTD 
I? (REX. LE. 0.0. AND.AIHX.LT. 0.) PX=1B0. +ATAN (AIHX/REX) *RTD 
IF (REX.GE. 0.0. AND.AIHX.LT. 0.) PX=360.-ATAN (AIHX/REX) ♦RTD 
PHASEX (H,N) =PX 
REY=REAL (Y(H,N)) 

AIRY=AIMAG(Y(M,N)) 

IF (REY .GE. 0.0. AND. AIRY. GE. 0.) PY=ATAN (AIHY/REY) *RTD 
IF (REY . LT. 0.0. AND. AIHY.G5. 0.) PY = 180.-ATAN (AIHY/REY) *RTD 
IF (FEY. LE. 0.0. AND.AIHY.LT. 0.) PY=180.*ATAN (AIM Y/RE Y) *RTD 
IF (REY. GE. 0.0. AND. AIHY.LT. 0. ) PY=360.-ATAN (AIHY/P.EY) *RTD 
PRASEY (H,N)=PY 
725 CONTINUE 

GO TO 900 

C **** NOW TAKE THE FOURIER TRANFORM OF X AND Y AND MULTIPLY BY Z 
C »** TO OBTAIN THE SCATTERED FIELDS * 

730 CALL FFT3D(X,IX,IX,IX,IX # 1,69,IWE,RWK,CKK) 

CALL EFT3D(Y,IX f TX,IX, IX , 1 , 6 P, INK, RHK, CU K) 

DO 760 fl = 1 , IX 
DO 750 N= 1 , IX 
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CXMN=X(N,N) 

X (PI , N) = ( (0 (M,N)-V(N,H) **2/G(M,N) ) *X(N,N) - ( O ( nj *V (N,N) /«{«, N)) 
.*Y (H,HJ)/(J*tf*2?)/2. 

Y(H,N) = (-U(H)*V{M,N)/G(N,H)*CXMN«>(G(H,N) -u (M) * *2/G (M ,N) ) 

750 CONTINUE 

760 CONTINUE 

C ***** CALCULATE THE REFLECTION COEPICTSNTS ****** 

CREFX=X ( 1, 1) /FLOAT (IX) **2 
CREFX=Y ( 1 , 1 ) /FLOAT (IX) **2 

CEEF= (CRSFX*SIN (-PHI/RTD) ♦CREFY*COS (PHI/RTO) ) 

CRET= (CREFX*COS (PHI/RTO) ♦CREFY*SIN (PHI/RTO) ) 

IF (ITil.GT.O) GO TO 800 
c *** TE POLARIZATION ***** 

C ** THIS IS THE CO-POLARIZED COfIPONENT **** 

RSFF=CABS (CREF/EF) 

C *** THIS IS THE CROSS— POLARIZED COMPONENT ****** 

PEFT=CABS (CRST/EP) 

WRITE (3 ,770) R EFF, REFT 
770 FORMAT (3X,2F10. 5) 

GO TO 900 
800 CONTINUE 

C *** TN POLARIZATION ***** 

C ** THIS IS THE CO-POLARIZED COMPONENT **** 

P.ETT=CABS (CPET/ET) 

C *** THIS IS THE CROSS-POLARIZED COMPONENT ****** 

RETF=CABS (CREF/ET) 

WRITE (3,770) RBTT,RETP 
WRITE (3 ,170) ITER 
STOP 
END 


O 


900 



139 


8.8 LISTING OF THE S.D.C.G. METHOD FOR THIN STRIPS WITH 
THE CROSS-SHAPED UNIT CELL 


C 

C*** MICHAEL DROZD-CIJRISTOS CURRENT FORMULATION ******* 
c *****CONJ. GRAD. METHOD FOR CURRENTS ON A CROSS ****** 

C**** MINIMIZATION IN THE RANGE **** 

COMPLEX CONE,CZERO,CXMN,CREFY 
COMPLEX G (32, 32J/1024* (0. 0, 0.0) / 

COMPLEX Y( 32, 32) /1 024* (0.0, 0.0) / 

COMPLEX X(32, 32) /1024*(0. 0,0.0)/ 

COMPLEX YO (32,32)/1024* (0. 0 , 0 . 0 )/ 

COMPLEX XU (32, 32) /1 024* (0. 0, 0.0) / 

COMPLEX RX (32,32) /1024* (0.0, 0.0) / 

COMPLEX RY (32, 32) /1 024* (0.0, 0.0) / 

COMPLEX J, HXI, HYI,CWK(32) 

COMPLEX DY (32,32) /1024* (0.0, 0.0) / 

COMPLEX DX (32, 32) /1 024 * ( 0. 0, 0. 0) / 

COMPLEX TX(32,32)/1024*(O. 0,0. 0) / 

COMPLEX TY (32, 32) /1024* (0.0, 0.0)/ 

REAL K,K2,RWK{342) 

DIMENSION IWK(342),RR(300),CH (300) ,X1 (1024) ,Y1(1024),Z1 (1024) , 
. Z2 (1024) , AMPJ[32) ,RINBEX (32) 

REAL U (32)/32*0.0/ 

REAL V (32,32) /1024*0.0/ 

C *** AA=SPACING BETWEEN VERTICAL WIRES' ***** 
c *** BB=THICKNESS OF VERTICAL WIRES ***** 

C *** CC=S?ACING BETWEEN VERTICAL WIRES ***** 

C *** DD=THICKNESS OP VERTICAL WIRES ***** 

READ (1 ,22) AA,BB,CC,DD,F 
22 FORMAT (5E10.. 4) 

P=2. 998E+8 

C *** IOPT«0 FOR A CROSS **** 

C ****IOPT=1 FOR A PARALLEL GRID **** 

IO PT* 1 

IF (IOPT.GT.0) CC=1«500E+15 
IF (IOPT. GT. 0) DD=0. 0 
WRITE ( 3, 33) AA,BB,CC,DD,F 

33 FORMAT (' 0* , ' AA»',F8.4,' THICK. OF VER. WIRE** , F8. 4, * CC= », 
.F8.4,' THICK. OF HOR. WIPE* 1 , F8. 4, * PREQ* »,E10.4) 

READ(1 ,22) PHI,THI,PSI 
WRITE (3,55) PHT,THI,PSI 

55 FORMAT (' O' , ' PHI=',F7.1,' THETA* ' , F7. 1 , ' PSI=*,F7.1) 

C *** READ SAMPLING NUMBER *** 

READ (1,66) IX 

C *** ITN=0 FOR TE POLARIZATION **** 

C *** ITM=1 FOR TM POLARIZATION **** 

C *** READ ITM ***** 

READ (1,66) ITM 

C *** READ NUMBER OF ITERATIONS *** 

READ { 1 , 6 6) NO I 
66 FORMAT (13) 

PI=3. 141593 
PI2=PI/2. 
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TPI=6. 23318*) 

OV= 2. 9979562*8 
UD=4.S-7*PI 
RTD=57. 29578 
EP=8. 954E-12 
ETA=SQRT (UU/EP) 

J=CHPLX(0. 0,1.0) 

ITEE= 1 

CONE=CMPLX (1 . 0,0.0) 

CZERO=CMPLX {0. 0, 0. 0) 

. H=TPI*F 
AL AHB=C7/F 
AA-AA/ALAHB 

C *** DETERMINE CONDUCTING REGIONS **** 

HDOX= (IX/2*1-BB*IX/(AA*2) ) 

HOPX=(IX/2*H-BB*IX/(AA*2) ) 

HDOY* (IX/2*1-DD*IX/(CC*2) ) 

H0PY=(IX/2*1+DD*IX/(CC*2) ) 

IF (IOPT. GT. 0) MDOI=IX*1 

K=TPI/ALANB 

K2»K**2 

STSPR=SIN (THI/RTD) ♦SIN (PHI/BTD) *K 
STCPK® SIN (THI/RTD) *COS (PHI/HTD) *K 
C?S=COS (PSI/RTD) /SIN (PSI/RTD) 

77 CONTINDE 

C ♦*♦ DEFINE THE FLOQIJET COEFFICIENTS •♦• 

DO 200 H=1,IX 
IF (H . GT. IX/2* 1 ) GO TO 125 
O (H) =TPI* (H- 1) /AA-STCPK 
GO TO 127 • 

125 O (H) s*TPI* (H-IX-1) /AA-STCPK 
127 CONTINDE 

DO 190 N=1,IX 

IF (H. GT. IX/2*1 • AND. N.GT. IX/2* 1) GO TO 160 
IF(N.GT.IX/2*1) GO TO 150 
IF (N .GT.IX/2*1) GO TO 140 

7(H,N)«TPI*(N-1) /CC-TPI* (H-1) /\ A*CP5-STSPK 
GO TO 170 

140 7 (N,N)=TPI*(N-IX-1) /CC-TPI* (H- 1) /A A*CPS-STSPK 

GO TO 170 

150 7 (H, N) =TPI* (N-1) /CC-TPI* (H-IX-1) /AA*CPS-STSPK 

GO TO 170 

160 7(N,N) =TPI* (N-IX-1) /CC-TPI* (H-IX-1) /AA*CPS-STSPK 

170 IF (K2.GE.U (H) **2*7 (H,N) **2) G (H, N) »-J*SQRT (K2- (U (N)**2*7 (H, N) **2 

.)) 

IF (K2.LT.0 (H) **2*7 (M, N) **2) G (H, N) =-SORT (D (H) **2*7 (M , N) **2-K2) 

. *CONE 

190 CONTINDE 
200 CONTINDE 

IF (TTM.GT. 0) GO TO 210 

C *** INCIDENT FIELDS FOR TS POLARIATION **** 

EXI=SIN (-PHI/RTD) 

EYI=COS (PHI/RTD) 

!1XI=C0 S (PH I/FT D) *C0S (THI/RTD) /ETA 
HYI=SIH (PHI/RTD) *COS (THI/RTD) /ETA 
GOTO 261 



C **** INCIDENT FIELDS FOR TK POLARIZATION **** 

710 HX I-CO 5 (PHI/RTD) *COS (THI/RTD) 

RY 1^5 1?) (-PHI/RTD) *COS (THI/RTD) 

HYI^SIN (PHI/RTD-PI2) /ETA 
IlXI=COS (PHI/RTD-PI2) /ETA 
261 CONTINUE 

C *** GIVE AN INITIAL GUESS ***** 

DO 310 9= 1 , IX 
DO 300 N=1,IX 
X(M,H)=CZERO 
Y(N,N) =CZER0 

IF (9. LE.BUPX. AND. 9. GE. HDOX) GO TO 270 
IF(N.LE.9UPY. AND.N.GE.BDOY) GO TO 270 
GO TO 280 

270 EN0R9=EN0F9*EXI*EXI*EYI*EYI 
280 XU(M,N)=X(9,N) 

YU (9,N) = Y (9,N) 

300 CONTINUE 
310 CONTINUE 

WRITE (3,320) IX, NOT, BDOX, 9UPX, HDOY , 9 UP Y, IT 9 
320 FORMAT (• 0* , * SA9P POINTS* *, 12, « # ITER=*,I3,* DOWN PNT X=*,I3, 

.* OP PNT X* * ,13, ' DOWN PNT I=*,I3,» OP PNT Y=',I3,* IT9= *,I3) 
C **** WORK ON INITIAL GUESS*** 

C *** TAKE THE POURIER TRANFORH OF THE INITIAL GUESS *•* 

CALL PFT3D {XU, IX,IX, IX, IX, 1,69, XWK,RWK,CWK) 

CALL FFT3D (YU, IX, IX, IX, IX, 1,69,IWK,RWK,CWK) 

C *♦* NULTIPLY INITIAL GUESS WITH THE HATRIX Z ***** 

DO 360 9=1, IX 
DO 350 N= 1 , IX 
CXBN=XU(9,N) 

XU (H, N) = { (G (H, N) -V (H,N) ** 2/G (H, H) ) *XU (H,N) - (U (9) * V(9,H) /G (H,N) ) 

. *YU(H,N))/(J*W*EP)/2. 

YU (9, N) = (-U (9) *V (9, N) /G(9,N) *CX9N* (G (9, N) -U (9) **2/G (9,N)) 

. *YU(9,N))/(J*W*2P)/2. 

350 CONTINUE 

360 CONTINUE 

C **• TAKE THE INVERSE 07 FOURIER TR\ HS70R9 ***** 

CALL 7FT 3D (XU, IX, IX, IX, IX, 1,-69,IWK,HWK,CWK) 

CALL FFT3D (YU, IX ,IX,IX ,IX, 1 ,-69 , IWK, RWK,CWK) 

C *** DEFINE THE RESIDUAL VECTORS RX AND RY ***** 

^ DO 450 9=1, IX 
DO 440 N=1, IX 
PX(M,N)=EXI*XU(9,N) 

HI (B, N) =EYI* YU (9 , N) 

IF (9.LE. 9UPX.AND.H.GE. HDOX} GO TO 400 
IF (N.LE. 9UPY. AND.N. GE. 9DOY) GO TO 400 
RX(M,N) =CZERO 
RY (B,N)=CZERO 

400 ER ROR= ERROR *RX (9,N) *C0NJG (RX (9,N) ) *RY (9, N) *CONJG (RY (9, N) ) 

DX (B,N)=RX(B,N) 

DY (B,N)=RY(B,N) 

440 CONTINUE 
450 CONTINUE 
C 

C *** FIND THE FOURIER TRANS OF TtlE RESIDUAL **** 

C 
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CALL FFT3D(DX. 17 . TX, 17 . TX, 1 , 69 , I WK , R WK ,CWK ) 

CALL FFT3D(DY,IX,IX,IX,IX, 1 , 69, I WK, R MK ,CHK) 

C *** CALCULATE THE DIRECTION VECTORS DX f- DY .**** 

DO 540 11=1, IX 
DO 530 ?J = 1,IX 
CXMN=DX (M,N) 

DX(M,H) = (CONJG {G (M, N) -V (M,N) **2/G (B, N) ) *DX(B,E) - 
. CONJG {V(N,N)*tJ (M) /G (B,N)) *DY (B, N) ) /CO NJG ( J*B*EP)/2. 

DY (M, H) = (CONJG (-V (M, N) *U(B) /G(B,N)) *CXBN* 

. COM JG (G (B,N) -0 (B) **2/G (B,N)) *DY (B,N)) /CONJG (J*W*EP)/2. 

530 CONTINUE 
540 CONTINUE 

C **** NOW FIND THE INV. FOURIER TRAS. ,OR THE DIBEC. PONCTIONS ** 

C 

CALL FFT3D (DX, IX, IX, IX, IX, 1 , -69 , 1 WK, RWK.CWK) 

CALL FFT3D(DY,IX,IX,IX,IX, 1,-69, IWX, RWK.CWK) 

C *** STORE DX AND DY IN TX AND TY ****** 

DO 560 B=1,TX 
DO 550 N=1,IX 

IF (B.LE.BUPX. AND.fl.GE.BDOX) GO TO 545 
IP (N.LE.MUPY. AND. N.GE.BDOY) GO TO 545 
DX (B, N)=CZERO 
DY (H, N) = CZERO 
545 TY(B,N)=DY(B,N) 

TX (B,N)=DX (B,N) 

P3=F3*CONJG (DX (H , N) ) *DX (B, N) *CONJG (DY (B, N) ) *DY (H, N) 

550 CONTINUE 
560 CONTINUE 

585 CALL FTT3D (TX,IX,IX,IX,IX,1,69,IWK,BWK,CWK) 

CALL PFT3D (TY, IX, IX, IX, IX, 1 ,69 ,IWK,HBK,CWK) 

C *** MULTIPLY TX AND TY BY THE HATBIX Z ***** 

do 610 a=i;ix 

DO 600 N=1 , IX 
CXHN=TX (B, N) 

TX (H, N) = ( (G (8 ,N) -V (H, N) **2/G (H, N) ) *TX (H,N) - (U (H) *V (H, N) /G (fl, N) ) 
. *TY (B,N))/(J*W*EP)/2. 

TY (H, N) = (-U (H) *V (B»N) /G (B,N) *CXHN* (G (B, R) -U (B) **2/G (B, H) ) 

. *TY(B,N) )/(J*W*EP)/2. 

600 CONTINUE 

610 CONTINUE 

CALL FFT3D (TX, IX , IX, IX , IX, 1,-69 ,IVK,RWR,CWK) 

CALL PFT3D (TY, IX, IX, IX, IX, 1 ,-6 9, IWK, RWR, CWK) 

FI =0.0 

DO 666 B=1 , IX 
DO 666 N= 1 ,IX 

IF (B.LE.BUPX. AND. B. GE. BDOX) GO TO 666 
IF (R.LE.BUPY. AND. N.GE. BDOYJ GO TO 666 
TX(fl,N)=CZERO 
TY (H , N) =CZERO 

666 F1 = F1 ♦CONJG (TX (B,N) ) *TX(fl,N) ♦CONJG (TY (B,N) ) *TY (S,N) 

ITER=ITER*1 

C *** EVALUATE THE PARAMETER AN ***** 

AN=F3/F1 

CH (ITER) = (ERROR) /(ENORM) 

C *** CALCULATE THE ERROR ****** 

ERROR= ERROR- (F3**2/F 1) 
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ERR = CH (TT5R) * 1 00 
PR (ITER) =FLOAT (ITER) 

WRITE (7, *) CH(ITSR) ,RR (ITER) 

C *** IMPROVE PREVIOUS ITERATE FOR X AND Y **** 
no 760 M=1,IX 
DO 750 N= 1 , IX 
X(M,M)=X(M,N)*AN*DX(M,N) 

Y (J1,N) = Y (M,N) ♦AN*DY (M, N) 

750 CONTINOE 
760 CONTINOE 

C **** UPDATE THE P.X AND R Y AND STORE THEM IN TX AND TY **** 

DO 843 M=1,IX 
DO 833 N=1,IX 
RX (ff,N) =RX (H,N)-AN*TX(M, H) 

RY(H,N) =RY (H, N) -AN*TY(M, N) 

TX(M,N)aRX(M,N) 

TY(H,N)=RY(H,N) 

833 CONTINOE 
843 CONTINOE 
C 

CALL FFT3D (TX, IX, IX, IX, IX, 1,69,IWK,RWK,CWK) 

CALL EFT 3D (TY, IX, IX, IX, IX, 1,69, I OK, RWK,CWK) 

C ***♦ flOLTIPLY TX C TY WITH THE COHJG. TRANSPOSE OP HATHIX 2 *** 

DO 863 Hal, IX 
DO 853 N=1,IX 
CXHN=TX(H,N) 

TX (H, N) = (CONJG (G (H, N) -V (H,N) **2/G (H, N) ) *TX (H, N) - 
. CONJG (V (H,H) *0 (H) /G (M,N) ) *TY (H, N) ) /CONJG { J*W*EP)/2. 

TY (M,N) = (CONJG (-V (H,N) *0 ( H) /G (H, N) ) *CXHN* 

, CONJG (G(H,N) -0(H) **2/G (H,H) ) *TY (H,N) ) /CONJG (J*W*EP)/2. 

853 CONTINOE . 

863 CONTINOE 

CALL PFT3D (TX, IX, IX, IX, IX, 1,-69, IWK, RWK,CWK) 

CALL PPT3D (TY,IX,IX,IX,IX,1,-69,IWK, RWK, CWK) 

P2=F3 

F3»0.0 

DO 900 H® 1 t ix 
DO 900 H 3 1 , IX 

IP (H.LE.HOPX.AND.H.GE. HDOX) GO TO 870 
IF (N.LE.HOPY.AND.N.GE. HDOY) GO TO 870 
TX (H, N) aCZERO 
TY (H,N) aCZERO 

870 F3aF3*CONJG (TX (H,N) ) *TX(H f N) >CONJG (TY (H, N) ) *TY (H,N) 

900 CONTINOE 

C *** DEFINE THE PARAMETER BN *** 

BN=F3/F2 

C *** CALCOLATE A NEW ESTIMATE FOR THE DIRECTION VECTORS DX Z DY *** 
DO 964 M=1 , IX 
DO 954 N=1,IX 
DX (H,N)=TX (H,N) ♦BN*DX (M, N) 

DY (M,N)=TY (R, N) ♦BN*DY (M, N) 

TX (H,N)=DX(H,N) 

TY (M , N) =DY ( M , N) 

954 CONTINOE 
964 CONTINUE 

IF (ERR. LT. 0.001) GO TO 1000 
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IF (ITER.LE. NOI) GO TO 585 
CALL GENPT (RR,CH,NOI,0) 

C GO TO 125 

970 DO 980 R=1,IX 
DO 980 N=1 , IX 
I=(R-1)*IX*N 
Z1 (I) *CABS (X (M ,N) ) 

Z2 (I) = CABS (X (R ,N) ) 

XI (I)* (FLOAT (fl-IX/2) 5) /IX*AA*1.040 

XI (I}= {PLOAT (N-IX/2) 5) /IX*AA*1.040 

WRITE (8, *) X1(I) ,X1 (I) ,Z1 (I) ,22 (I) 

980 CONTINOE 
1000 DO 1100 1=1, IX 

ARP(I) =CABS (X (17,1) ) 

RIMDEX (I) = (PLOAT (I-IX/2) 5) /IX*AA*1.045 
1100 WRITE (8,*) ARP (I) , RIMDEX (I) 

STOP 

END 
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8.9 LISTING OF THE S.D.C.G. METHOD FOR APERTURE FIELDS 


C*****C0N.T. (jit AD. METHOD 2 .FORT***** 

C *** SOLVES FOE THE APERTURE FIELDS ******* 

C **** MINIMIZATION IN THE RANGE ********* 

COMPLEX CONE,CZERO,CXMN,CREFX,CREFY 
COMPLEX G (32, 32)/1024* (0.0, 0.0) / 

COMPLEX Y(32,32)/1024* (0.0, 0.0) / 

COMPLEX X (32,32) /I 024* (0.0, 0.0) / 

COMPLEX YU (32,32)/1024* (0. 0, 0.0)/ 

COMPLEX XO(32, 32) /1 024 *(0.0, 0.0) / 

COMPLEX RX (32, 32)/1024*(0. 0,0.0) / 

COMPLEX RY{32,32)/1024*(0. 0, 0.0) / 

COMPLEX J , HXI, HYI , CW K ( 32) 

COMPLEX DY (32,32)/1024*(0. 0,0.0) / 

COMPLEX DX(32,32)/1024*(0. 0,0.0) / 

COMPLEX TX{32,32)/1024*(0. 0,0.0) / 

COMPLEX TY (32, 32) /1 024* (0.0, 0.0) / 

REAL K ,K2, RWK ( 342) 

DIMENSION INK ( 342) ,RR(250) ,CH(250) ,Z1 (1024) ,X1 (1024) ,11 (1024) , 
.22(1024) , AMP (32) ,RINDEX (32) , CROSS (32) ,PHASEX (32,32) , PHASE! (32,32) 
REAL U (32) /32*0. 0/ 

REAL V (32,32) /1 024* 0.0/ 

INTEGER COUNT 

HEAD (1,22) AA, BB,CC, DD, F, ERR 
22 FORMAT (8S10. 4) 

F=2.998E*8 

C *** IOPT=0 FOR A SQUARE OR A RECTANGULAR MESH ****** 

C *** IOPT*1 FOR A PARALLEL GRID ****** 

IOPT*0 

IF (IOPT.GT. 0) CC*1. 500E+15 
IF (IOPT.GT.0) DD=1 . 500E* 15 
HRITE(3,33) AA, BB, CC ,DD, ERR 
33 FO RM AT ( • 0 * , • A* ',F15.8,* B= *,F15.8,» C= *,F15.8, 

9* D= *,P15.8,» ERR* »,F15.8) 

WRITE (3, 44) F 

44 FORMAT (*0* ,’ FREQ * *,E10.4) 

READ (1 , 22) PHI,THI,PSI 
WRITE (3, 55) PHI, THI, PSI 

55 FORMAT (' 0* , • PHI* ',F10.1,» THETA* ',F10.1,» PSI* »,F10.1) 

C *** READ NUMBER OF SAMPLING POINTS *** 

READ (1,66) IX 
READ (1,66) ITM 
READ (1,66) NOI 
WRITE (3,56) ITM 

56 FORMAT (3X, 'THE VALUE FOR ITM IS* ',13) 

66 FORM AT (13) 

PI=3. 1 41 593 
PI2=PI/2. 

TPI*6. 283185 
CV=2.997956E+8 
Un=4.E-7*PI 
P.TD*57. 29578 
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E? = 8.8 54E-12 
ETA=5QRT (00/ EP) 

J=CMPLX(0.0, 1.0) 

ITEP = 1 

cone=crpix (i.o ,o.O) 

CZERO=C»PLX (0. 0, 0. 0) 

W=TPI*F 

ALAHB=CV/F 

AA=AA/ALANB 

BB=BD/ALAHB 

CC=CC/ALANB 

DD-DD/AL ABB 

NX=IFIX(BB/AA*PLOAT(IX) *2.) /4*2 

NT =IPIX (DD/CC* FLOAT (IX) *2. ) /4*2 

NX 1= (IX- NX) /2+1 

NX2=NX1*NX-1 

NT 1= (I X-NY) /2* 1 

NY2=NYUNY-1 

WRITE (3, 70) NX,NX1,NX2 # NY,NY1,NY2 
70 FOBS AT (' 0* , • NX=*,I3,3X, , NX1=»,I3,3X,*NX2=»* ,13, 3X, 

. • NY=* ,I3,3X,*NY1=* ,I3,3X,'NY2=* ,13) 

K=TPI/ALAH8 

K2=K»*2 

STSPK=SIN (THI/RTD) *SIN (PHI/BTD) *R 
STCPR-SIN (THI/RTD) *COS (PHI/RTD ) *K 
CPS=COS (PSI/HTD) /SIN (PSI/RTD) 

75 CONTINOE 

c «** deterhine floqoet coefficients **** 

DO 120 H=1,IX 
IF (R.GT. IX/2*1) GO TO 80 
0 (N) =TPI* (H— 1) /AA-STCPK 
GO TO 85 • 

80 0 (H) ®TPI* (H-IX-1) /AA-STCPK 

85 CONTINOE 

DO 115 N=1 , IX 

IF (H.GT. IX/2+1 . AND.N.GT. IX/2+ 1) GO TO 100 
IP (H.GT. IX/2+1) GO TO 95 
I? (N.GT. IX/2+1) GO TO 90 

V(H,N) =TPI*(N-1) /CC-TPI*(H-1)/\ A*CPS-STSPK 
GO TO 110 

90 V (H,N)=TPI* (N-IX-1) /CC-TPI*(H-1) /AA*CPS-STSPK 
GO TO 110 

95 V(R,N) =T?I* (N-1) /CC-TPI* (H-IX-1) /A A*CPS-STSPK 
GO TO 110 

100 V(H,N) =TPI* (N-IX-1) /CC-TPI* { H-IX-1 ) /A A*CPS-STSPK 
110 IF (K2. GE.O { N) **2*7{H,N) **2) G(H,N) =-J»SQRT (K2- (0 (H)**2*V (fl,N)**2 
.)) 

IF (K2. LT.O (M) **2*V (N,N) **2) G (M, N) =-SQRT (0 (R) **2*T (H,N) **2-K2) 

. *CONE 

115 CONTINOE 
120 CONTINOE 

IF(ITH.GT.O) GO TO 130 

C *** INCIDENT FIELDS FOR TE POLATIZATION **** 

EXI= SI N (-PHI/RTD) 

ETI=COS (PHI/RTD) 

HXI=COS (PHI/RTD) *COS (THI/PTD) /ETA 
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HY T= 5111 (PHI/RTP) *C05 (TliT/B TD) /ETA 

r, ore mo 

C *** INCIDENT FIELDS FOB TM POLARIZATION *** 

130 F.XI=COS (PHI/P.TD) ’•COS (THI/BTD) 

EYI=SIN (-PHI/RTD) *COS (THI/RTD) 

HYI=SIN (PHI/RTD-PI2) /ETA 
HXI=COS (PHI/RTD-PI2) /ETA 
140 CONTINUE 

C *** GIVE AN INITIAL GUESS ****** 

DO 145 M=NX 1 , NX2 
DO 142 N=HY1,NY2 
X (M,N) =CZERO 
Y(M,N) =CZEHO 
XU (M,N)=X (M,N) 

YU (M,N)=Y(M,N) 

F6=F6*CONJG (HXI) *HXI*CON JG (HYI) *HYI 
142 CONTINUE 
145 CONTINUE 

C **** WORK ON INITIAL GUESS*** 

C **** MULTIPLY INITIAL VECTORS XU 6 YU BY THE MATRIX Z *** 

C 

C *** TAKE THE POURIER TRANSFORM OF XU B YU ***** 

CALL FFT3D(XO,IX,IX,IX,IX, 1,69,1 WK,RWK,CWK) 

CALL FFT3D(¥U,IX,IX,IX,IX, 1,69,IIK,RWK,CWK) 

DO 160 H*1,IX 
DO 150 N=1,IX 
CXMN=XU(M,N) 

XU(H # N) = (U(M) *V (H,N)/G(N, N) *XU(N,N) ♦ (V (M, N) **2/G{H,N)-G (H,N) » 
. *YU (H,N) ) * (J/W/UU) 

YU (M, N) = { (G (M ,N) -U (M) **2/G (M, N) ) *CXHH-U (H) *V (H, N) /G (H, N) 

. *YU (H,N) ) * (J/W/UU) 

150 CONTINUE • 

160 CONTINUE 

CALL FFT3D(XO,IX,IX, IX, IX, 1,-69,IWK,P.WK,CVK) 

CALL FPT3D(YU,IX,IX,IX,IX,1,-69,IWK,BWK,CWK) 

C *** CALCULATE THE RESIDUAL VECTORS RX 6 RY *** 

DO 190 H=1,IX 
DO 180 H»1,IX 
RX (H,N)=HXI-XU (H,N) 

RY (M,N)=HYI-¥U (M,N) 

IF (M. GE. NX 1 « AND. M. LE. NX2. AND. N.GE. NY 1. AND. N. LE. HY2) GO TO 175 
RX (M,N) =CZERO 
RY (M,N)=CZERO 

175 ERROR=ERROF*RX {M f N) *CONJG (RX (M, N) ) ♦RY (M, N) *CONJG (RY (N, N) ) 

DX (M,N)=RX(M,N) 

DY (N,N)=RY(M,N) 

180 CONTINUE 
190 CONTINUE 

**** MULTIPLY THE RESIDUALS BY THE CONJG. TRANS. OF Z 

TO FIND THE DIRECTION VECTORS DX & DY ***** 

**** FIND THE FOURIER TRANSFORM OF RESIDUALS *** 

CALL FFT3D(DX,IX,IX,IX,IX, 1 , 69 ,1 WK ,RWK,CVK) 

CALL FFT3D (DY, IX, IX, IX, IX, 1,69,IWK ,RVK ,CWK) 

DO 210 M=1 , IX 
DO 200 N = 1 , I X 
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CXBN=DX(B,N) 

DX (M, N) = (CONJG { 0 ( «) *V (S,N)/G (M, N) ) *DX (B,N) ♦ 

. CONJG { (G (B,N) -U («) **2/G (B,N) ) ) *DY (H,N) ) *CONJG (J/W/UU) 

DY (B, N) = (CONJG { (V (B,N) **2/G (B, N) -G(B,N) )) *CXBN- 
. CONJG (U (B) *V (B , N) /G (B, N) ) *DY (B, N) ) *CONJG (J/W/UU) 

200 CONTINUE 
210 CONTINUE 

C **** NOW FIND THE INV. FOURIER TEAS. ,OR THE DIREC. PONCTIONS ** 

C 

CALL FFT30 (DX, IX, IX, IX, IX, 1,-69, IWK, RWK, CWK) 

CALL' PPT3D(DY,IX,IX,IX,IX,1,-69,IWK,RWK,CWK) 

C *** STORE DX & DY IN TX 6 TY ****** 

DO 230 B= 1 , IX 
DO 220 N=1,IX 

IP (B.GE. NX 1. AND. B. LE. NX2 .AND. N.GE.NY1.AND. N. LS. NY2) GOTO 215 
DX (B,N)=CZERO 
DY (B,H) = CZERO 
215 TY(B,N)=DY(B,N) 

TX (B,N) = DX (B,N) 

F3=F3*C0NJG(DX (B ,N) ) *DX(B,N) +CONJG (DY (B,N) )*DY (B,N) 

220 CONTINUE 
230 CONTINUE 

C *** THE ITERATIVE PROCESS STARTS NOW S I ! t •*** 

C **** BULTIPLY THE DIRECTION VECTORS BY THE BATRIX Z ***** 

240 CALL PFT3D (TX, IX, IX, IX, IX, 1,69,IWK,RWK,CWK) 

CALL FFT3D (TY, IX,IX, IX,IX, 1 , 69,1 WK , RWK , CWK) 

. DO 261 B=1,IX 
DO 251 N=1 , IX 
CXBN®TX(B,N) 

TX (B, N) = (U (B) *V (B ,N) /G (B , N) *TX (B, N) ♦ (V (B, N) ** 2/G ( B, N) -G (B,N) ) 
. *TY (B ,N) )'* J/W/UU 

TY (B, N) = ( (G (B,N)-U (B) **2/G(B, N) ) *CXBN-U (B) *V (B,N) /G (fl, N) 

. *TY(B,N))* J/W/UU 
251 CONTINUE 

261 CONTINUE 

CALL FFT3D (TX, IX, IX, IX, IX, 1,-69, IWK, RWK, CWK) 

CALL FFT3D (TY, IX, IX, IX, IX, 1, -69, IWK , RWK, CWK) 

P1=0.0 

DO 300 B= 1 , TX 
DO 300 N=1 , IX 

IF(B.GE.NX1.AND.B.LE.NX2.AND.N.GE.NY1. AN D. N. LE. NY2) GO TO 300 
TX (B,N)=CZERO 
TY (B,N) = CZERO 

300 F1=F 1+CONJG (TX(H,N)) *TX (B, N) *CONJG (TY (B, N) ) *TY (B,N) 
ITEF.=ITER*1 

C **** CALCULATE THE FACTOR AN **** 

AN=F3/F 1 

CH(ITER) =SQRT( ERROR) /SQRT ( F6) 

C **** CALCULATE THE ERROR **** 

ERROR=ERROR- (F3**2/F1) 

C *** UPDATE THE VALUES FOR X f> Y ***** 

DO 410 B=1 , IX 
DO 400 5*1, IX 
X (B, N) = X ( B , N ) *AN*DX (B, N) 

Y (B, H) =Y (M , N) ♦AN*DY{B,N) 

400 CONTINUE 
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a io continue 
F5=0. 0 

C *** FIND A NEW ESTIMATE FOR THE RESIDUAL VECTORS RX S RY *** 

DO 443 M=1,IX 
DO 433 N=1,IX 
RX (fl # N)=RX(M,N)-AN*TX(M,N) 

HY (3,N)=BY (N,N) -AN*TY(M,N) 

TX (8,N) =RX (H,N) 

TY(H,N)=RY(S,N) 

433 CONTINUE 
443 CONTINUE 

RF. (ITER) = FLOAT (ITER) 

C WRITE (8 ,*) CH (ITER) , HR (ITER) 

C *** HUTLTIPLY TX & TY BY THE CONJG. TRANS. OP THE MATRIX Z ♦** 
CALL PFT3D (T X, IX ,IX , IX , IX, 1 , 69,1 WK, RWK # CHK) 

CALL PFT3D (TY,IX,IX,IX,IX, 1 , 69 ,IWK, HWK,CWK) 

DO 460 H=1,IX 
DO 450 N=1,IX 
CXBN*TX (H,N) 

TX (H,N) = (CONJG (U(B) *V(M,N)/G(H, N) ) *TX (H,N) ♦ 

CONJG ( (G (8, N) -U (N) **2/G (S,N) ) ) *TY (H,N) ) *CONJG (J/W/UO) 
TY (f!, N) = (CONJG ( (V (H,N) **2/G (H, N) -G (H, N) ) ) *CXHN- 
. CONJG (U(H) * V ( H , N) /G (8, N)) *TY (H,N) ) *CO NJG ( J/W/UU) 

450 CONTINUE 
460 CONTINUE 

CALL PFT3D (TX, IX ,IX, IX, IX , 1,-69,IWK # RWK f CWK) 

CALL FFT3D (TY,IX, IX, IX, IX, 1,-69, IWK,RWK,CRK) 

F2=F3 
F3=0. 0 

DO 470 8=1, IX 
DO 470 N= 1 , IX 

IF (8.GE. NX1 . AND. 8.LE.NX2.AND. N. GE. NY1. AND. N. LE.NY2) GO TO 465 
TX (H,N)=CZERO 
TY(H,N)=CZERO 

465 F3=F3 ♦CONJG (TX (8, N) ) *TX (8 , N) *CONJG (TY (8,N) ) *TY(B,N) 

470 CONTINUE 

C *** CALCULATE THE FACTOR BN ***** 

BN=F3/F2 

C **** UPDATE THE DIRECTION VECTORS DX S DY **** 

DO 564 8*1, IX 
DO 554 H=1 , IX 
DX (H,N) =TX (B,N) *3N*DX(8, N) 

DY (M,N) =TY (M,N) ♦ BN*DY(B, N) 

TX (N,W)=DX(3,JI) 

TY (N , N) = DY (3,N) 

554 CONTINUE 
564 CONTINUE 

C *** GO FOR ANOTHER ITERATION IF YOU WANT **** 

IF (ITER. GT. NO I) CALL GENPT (RR,CH,NOI,0) 

IF (ITER. GT. HOI) 00 TO 800 
IF (ERROR.LT. 0.0001) GO TO 800 
GO TO 240 

C *** STORE X S Y INTO THE 1-DIH. ARRAYS Z1 C Z2 TO BE USED POR 
C AMY PLOTTING PURPOSES 

570 DO 500 3 = 1 , IX 
DO 50 0 N = 1 , I X 


** * 
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1= (H-1)*IX*N 

7.1 (I)=CABS{X(H,N)) 

Z2 (I) = CABS (Y (H,N) ) 

XI (I) = (FLOAT (H-IX/2) -. 5) /IX*AA*1. 05 
Y1 (I) = (FLOAT (N-IX/2) -.5) /IX*BB*1.05 
WRITE(7,*) XI (I) ,Y1 (I) ,Z1 (I) ,Z2 (I) 

590 CONTINUE 
500 DO 725 »!=NX1,NX2 

DO 725 N=NY1,NY2 
BEX=REAL(X(H,N)) 

AIHX=AIH AG (X (H ,N) ) 

IF (REX .GE. 0.0. AND. AIHX.GE. 0.) PX=ATAN (AISX/RBX) *BTD 
IF(REX.LT.0.0. AND. AIHX.GE. 0.) PX=180.-ATAN (AIHX/REX) *RTD 
IF (REX.LE. 0.0. AND. AIMX. LT. 0. ) PX=180.*ATAN (AIHX/BEX) *HTD 
IP (R EX. GE. 0. 0. AND. AIHX.LT. 0. ) PX=360.-ATAN (AIHX/REX) *BTD 
PHASEX (H,N) =PX 
BEY=REAL(Y(H,N)) 

AIHY=AIHAG (Y (H,N) ) 

IF (BEY.GB. 0. 0. AND. AIHY.GE. 0.) PY=ATAN (AIHY/REY) *BT D 
IF (REY.LT. 0. 0. AND. AIHY. GE. 0. ) PY=1 80 . -ATAN (AIHY/REY) *RTD 
IF (BEY. LE. 0.0. AND.AIHY.LT. 0.) PY=180. + ATAN (AIHY/REY) *RTD 
IF (REY.GE. 0.0. AND.AIHY.LT. 0.) PY=360.-ATAN (AIHY/BE Y) *BTD 
PHASEY (H,N)=PY 
725 CONTINUE 
800 DO 820 1=1, IX 

AHP(I) =CABS (Y ( I, 16) ) 

R INDEX (I) = (FLOAT (I-IX/2) -.5) /IX*AA*1.0 45 
WRITE (8,*) AHP(I) ,SINDEX (I) ...... ... 

C CROSS (I) =CABS (Y (9, 1) ). 

820 CONTINUE 

CALL GENPT (RINDEX ,AHP, IX, 0) 

C CALL GENPT (RINDEX, CROSS, IX, 0) 

WRITE (3,840) ITER 
84 0 FOP.H AT (3 X, 13) 

900 STOP 
END 


nnnnnonnnnnnoonnnoonnnnnnnnnnnoononnnoonnnoooononnnoonnoonnnononon 


8.10 


LISTING OF ONE, TWO, AND ‘THREE DIMENSIONAL COMPLEX 
FAST FOURIER TRANSFORM 


1251 RCUTI5Z NAME 


FFTCC 


- I5M/B0UBLE 

- JANUARY 1, 1973 

- COMPUTE THE EAST FOURIER TRANSFORM OP A 
COUPLE! VALUED SEQUENCE 

- CALL PPTCC (A,H,IWK,WK) 

- COMPLEX VECTOR OP LENGTH N. ON INPUT A 
CONTAINS THE COHPLEX VALUED SEQUENCE TO BE 
TRANSFORMED. ON OUTPUT A IS REPLACED BY THE 
POORIER TBANSFORS. 

- INPUT NUMBER OP DATA POINTS TO BE 
TRANSFORBED. H BAY BE ANY POSITIVE 
INTEGER. 

- INTEGER WORK VBCTOH OP LENGTH 6*N*150. 

(SEE PROGRAMMING NOTES FOB PURTHEB DETAILS) 

- REAL NORN VECTOR OP LENGT! 6*N*150. 

(SEE PROGRAHNING NOTES FOR PURTHER DETAILS) 

- SINGLE AND DOUBLE/H32 

- SINGLE/H36, H48, H60 

8EQD. INSL ROUTINES - NONE B SQUIRED 

NOTATION - INPORNATION ON SPECIAL NOTATION AND 

CONVENTIONS IS AVAILABLE IN THE NANUAL 
INTRODUCTION OR THROUGH INSL ROUTINE UHELP 

BEN ABES 1. PPTCC COBPUTES THE FOUBIER TRANSFO* H, X, ACCORDING 
TO THE FOLLOWING FORMULA; 

X (K* 1) = SUB PBOH J * 0 TO N-1 OF 

A (J*1) *CEXP ( (0. 0. (2.0*PI«M*K) /N) ) 

FOR K=0, 1 , , N- 1 AND PI=3. 1415... 

NOTE THAT X OVERWRITES A OR OUTPUT. 

2. PPTCC CAN BE USED TO COMPUTE 


COMPUTER. 

LATEST REVISION 
PURPOSE 

USAGE 

APS UN ENTS A 

S 

IV FC 
NX 

PPECISIOH/HAEDHAPE 


X(K*1) 


0 TO N-1 OF 
) *CEXP ( (0. 0. (-2. 0*PI*J*K) /N) ) 


(1/N)*SUR PBOB J 
A (J + 1) *CEXP ( (0. 0. K ~*. 

FOR 1C=0, 1 , N— 1 AND PI=3.1415... 

BY PERFORMING THE FOLLOWING STEPS; 

DO 10 1=1, N 

CON JG (A (I) ) 

10 CONTINUE 

CALI. FFTCC (A , N , IWN, WR) 

DO 20 1=1, n' 

A inlir* COf,JG < ft <i>>/N 


COPYRIGHT 

VARPANTT 


20 CONTINUE 

- 1978 BY I HSL, INC 


ALL BIGHTS RESERVED. 


INSL WARRANTS ONLY THAT INSL TESTING HAS BEEN 
APPLIED TO THIS CODE. NO OTHER WARRANTY, 
EXPRESSED OR IMPLIED, IS APPLICABLE. 


SUBROUTINE FFTCC (A,N,IWK,WK) 

SPFCIFTC»TTONS FOR * BGUHENTS 



ORfGfMAl. VkQ£- *S 

Of POOR QUALITY 


INTEGER 

DOUBLE PRECISION 
CONPLEX*16 

INTEGER 


1 

2 

3 

1 

1 

2 


DOUBLE PRECISION 


COMPLEX*! 6 
EQUIVALENCE 


10 


15 


20 


DATA 

I 

DATA 


I? 
K = 
H = 
J a 
33 
3 ? 


(N 

N 

0 

2 


.EQ. 1) GO TO 9005 


= 1 . 


N.INK ( 1) 

HK Cl) 

A (N) 

SPECIFICATIONS POR LOCAL VARIABLES 
T.IAN.IAP.IBM. IBP. IC, ICC.ICF , ICK . ID, IDM1 , II. 
IJA. IkB, IRT, ILL, lft,IRD,ISP.ISK,ISP,ISS,ITA,ITB, 
J,JA,JF,JJ. JK,K,K0,K1,K2,K3,RA,KB,KD2,KF,KH,Kir; 
KT,Kfp.L,L1 # s!nft.E&1,ftp 

CM,SN,C1,C2,C3 # S1,S2,S3,C30,RAD,A0,A1.A4.B9, 
A2,A3,BO,B1,B2,B3,ZESO,HALF,ONE,TNO,ZO (2) , 

ZAO-ZAI.ZAZ^ZA j.2A4 # AK2 
5*0,Z0tl)) # (ZA^ZI^I))., 

Z Aj 
B 1 

B3, _ 

ftAD/6 

C30/.86fi02540378443S6DO/ 

ZER0,HAL7, ONE, THO/O. 0D0, 0. 5D0, 1. 0D0.2. 000/ 
FIRST EXECUTABLE STATE SENT 



INK (1) 

I =* K/JJ 
IP (I*JJ .NE. K) GO TO 10 
S a H*1 
INK (S+ 1) a 3 
K a I 
GO TO 5 
3 = 3*2 

IP (J .SQ. 9) J 3 3 

JJ a j * j 

IP (JJ .LE. K) GO TO 5 
KT = H 


DETER SINE THE SQUARE FACTORS OP N 


3 = 2 
I 3 K / 

IF <I*J 
S a H ♦ 

IRK (H*1) a 
K a I 

GO TO 15 
3 » J ♦ 1 
IP (J .EO. 
J a J ♦ 1 
IP (J.LB.K) 


DETERHINE TBE REMAINING FACTORS OP I 


25 

30 


K a INKfS+1) 

IP (IWKjKT* t 
IP (KT.LE.O) 

KTP a KT ♦ 

DO 25 I = 

J a KTP 

S a H*1 

IWKfH* 1) 
CONTINUE 
HP = H+1 
IC = HP*1 
ID = IC*NP 
ILL = I D* HP 
IRD = ILL*N?+1 
ICC = IRD+SP 
ISS = ICCtNP 
ICK = I55*flP 


NE. K) GO TO 20 


3) GO TO 15 
GO TO 15 


) .GT. INK (S+ 1) ) 
GO 


TO 30 

2 

1, KT 

= INK (.1) 


INK (KT*1 ) 
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ISK 

= 

ICK*K 


ICF 

S 

I SK *K 


ISF 

a 

ICF*K 


TAP 

s 

ISF+K 


KD2 

— 

(K“ 1) 

/ 2 

TBP 

a 

I AP ♦ 

KD2 

TAH 

a 

IBP ♦ 

KD2 

IBM 

a 

I AM ♦ 

KD2 


35 


40 


IWK ( J- 1) ♦ 
UK(ILL*L) . 


= 0 


5J«fl 


IWK 

2Q 


) .eq. 

GO TO 


♦ 1 


HM1 = M-1 
I = 1 

L = HP - I 

j = ic - r 

IWK(ILL*LJ_ = 0 

IP 

IP , 

1 = 1*1 
1 = L - 1 
IWK (ILL*L) 

I = I * 1 
IP(I.LE.HHl) GO TO 35 
— a 'ILL*1) = 0 
ILL* HP) = 0 
iq = i 
ID) = N 
5 J = 1,H 
K = IWK { J* 1) 

IHK(IC*Jl = IWK (IC*J-1) * 
IWK ?ID*J) =» IWK(ID*J-1) / 
WK (IRD*J) = 5AD/IWK (IC+J) 
Cl = RAD/K 

LE. 21 GO TO 45 
faCOS(CI) 


0 4 




IWK (ILL*L) = 1 


IWK 

IWK 

IWK 

IWK 

DO 


IP (K 
(ICC*J) 


WK 

nc = Dsiw(ci) 

45 CONTINOE 
HH 
IP 
IP 
sn 

CH 
- SH 
50 KB 
KH 


fIWK(ILL*H) . EQ. 
(HH .LB. ij GO TO 
= IWK (IC*HH-2) 

= DCOS (SM) 

= DSIM(SH) 

= 0 
= N 
0 


V" 

WK(IRD*f») 


= n - i 


55 


JJ 

I = 1 
Cl = OWE 
SI = ZEPO 
LI = 1 

IF (IWK (II.L*I*1) 
KF = IWK (1*1) 


, EQ. 1) GO TO 60 


C 

C 


GO TO 65 
60 KP = 4 
I * 1*1 

65 ISP = IWK (ID* I) 

IF (LI . EQ. 11 GO 
SI = JJ * WK (I ~ 

Cl = DCOS 

SI = DSIN 


m 


TO 70 
ED* I) 


JO 

75 KO 
K2 
IF 
80 KO 
IF 
K2 


IP (KF .GT. 41 GO TO 140 
GO TO (75,75,90,115) , KF 
= KB ♦ T S? 

= KO ♦ ISP 
(LI . EQ. 1) GO TO 85 
= KO - 1 

(KO .LT. KB) GO TO 190 
= K2 - 1 


FACTORS OF 2, 3, AND 
HANDLED SEPARATELY. 


ZA4 s A(R2t]) 


ARE 
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85 


90 


95 


AO a A4*C1-B4*S1 
BO = A4*51*B4*C1 
A (K2*1> = A (KO* 1) -ZAO 
AjK0*ll = A (KO* 1) *ZA0 
TO 80 


GO 

KO 

IP 


KO - 1 


(KO .LT. KB) GO TO 190 
K2 = K2 - 1 
AK2 a A (K2*1) 

K2*1) = A (K0* 1) -AK2 
KO* 1 j = A (KO* 1) *AK2 
TO 85 

IP (LI . EQ. 1) GO TO 95 
C2 = Cl * Cl - SI * SI 
S2 * TWO * Cl • SI 
JA = KB ♦ ISP - 1 
KA a JA ♦ KB 
I KB a KB*1 


an 

U 


IJA = JA* 1 
DO 110 II = 


IKRjI JA 


100 

105 


110 

115 


120 


K1 » k8 ♦ I5p 

K2 a K1 ♦ ISP 
ZAO = A (K0*1) 

IP (LI .EQ. 1) GO TO 100 
ZA4 a A(KT*1) 

A 1 = A4*C1-B4*S1 
B1 a A4*S1*B4*C1 
ZA4 a A (K2*1) 

A2 a A4*C2-B4*S2 
B2 a A4*S2*B4*C2 
GO TO 105 
Z A 1 a A(K1*1] 

ZA2 


> fK 1* 1 ) 
(K2*l[ 
a DC MI 


A(K0*1) = DC 
AO = -HAL? * 
A 1 
BO 

B1 a ( 

A (K1*1 
A (K2*1 
1IBE 
190 

(Li .EQ 


CONT 
GO TO 


iPLX (A0*A1*A2,B0*B1*B2) 
(A1*A2) ♦ AO 
( A 1-A2) * C30 
-HALF * (B1+B2) ♦ BO 
B1-B2) * C30 
) = DCMPLX (AO-B1,BO*A1) 

j = dcmplx (ao*bi ,bo-ai 


IP (L1 .EQ. 1) GO TO 120 
C2 a Cl * Cl - SI * SI 

52 a TWO * Cl * SI 

C3 a Cl * C2 - SI * S2 

53 = SI * C2 ♦ Cl * S2 
JA = KB ♦ ISP - 1 

KA = JA ♦ KB 
I KB = KB*1 
IJA = JA* 1 
DO 135 II = IKB,IJA 
KO = KA - II ♦ 1 
K1 = KO * TSP 
K2 a K1 ♦ TSP 
K3 a— Z 2 ♦ ISP 
ZAO = A (K0*1) 

IP (LI .EQ. 1) GO TO 125 
ZA4 = A (K 1 ♦ 1) 

A 1 = A4*C1-B4*S1 
B 1 = A 4*S 1 *B 4*C1 
7 , A 4 a ft (K2*1) 

A2 = A4*C?-B4*S2 
02 = A4*S2*BU*C2 
7A4 a A (K3* 1 1 
A3 = A4*C3-B4*S3 
B1 a au*?j3*RU*c3 



155 


125 


130 


135 

140 


145 

150 


155 

160 


GO TO 130 
ZA1 = A (K1+1J 
ZA2 = A/K2+1) 
Z A3 a A (K3* 1 i 
A(K0*1) 


A (K1 *1. 

A (K2*1j 
A]k3*1 

cohtimoe 

GO TO 190 
= FTP - 
= JK/2 
* ISK(IO *1 
= KB ♦ IS? 


= DC?1PLX (AO* \2*A1*A3,B0+B2*31*B3) 
= DCNPLX A0*A2-A1-A3 # B0*B2-B1-B3 
= DCflPLX A0-A2-Bl*B3,n0-B2*Al-A3 
a DCKPLX AO-A2*B1-B3,BO-B2-Al*A3) 


1 


JR 
KH 

K3 * iwk7id*i~i ) 

KO = KB ♦ 

IF (LI .SO 
K * JK - 1 
WK (ICF* 1) => Cl 
UK (ISP*1) a SI 
DO 145 J = 1 - fC 
WKIICP*J*1 
WK ?l SF ♦ J*1 
CONTINUE 
IF (KF .EQ. JF> 
C2 = WKfICC*I1 
WK(ICK*1) = C2 
WK (ICK* JK) = C2 
S2 = WK(ISS*I) 
WK(ISK*11 = S2 
WK(lSK*JK) * -S2 
DO 155 J - 1.KH 
K = JK - J 
WK(ICK*K1 
WK TCK*J*1 
WK ISK*J*1 
iK ISK+K) 
CONTINUE 


1) GO TO 150 


1 = WK 

[ICF*J1 

* Cl 

- WK 

[isf*ji * 

) a fc’K 

(ICF*J) 

* SI 

♦ WK 

[ISF*J) * 


SI 

Cl 


GO TO 160 


WK (ICK*J) * C2 
a wK (ICK*K) 


- WK (ISK*J) ♦ 

a WK(lCK*jJ. * S2 ♦ WK(ISK*J) * 


S2 


= -WK(lSK*J*1) 


C2 


KO = 
K 1 
K 2 
ZAO 


A3 

B3 

DO 


KO 

KO 

KO 

A 


- 1 

♦ K3 
(KO* 1) 


= AO 

17^' J a 1 KH 

K 1 a K 1 ♦ ISP 
K2 a K2 - ISP 
IF (11 .EQ. 1) 
KF * J 


GO TO 165 


ZA4 = A (K1 *1 ) 

A 1 a A4*WK(TCF*J)-B4*WK (ISF+J) 
B 1 a A4*WK (TSF*J) *B4*WK (ICP*Jj 
ZA 4 a A (K2*1 ) 

A2 a A4* WK (ICF+K) -B4*WK (ISF*K) 
B2 = A4*WK(ISF*K) *B4*WK (ICF*K) 
„ GO TO 170 

165 ZA1 = A (K1 *1) 

ZA2 a A(R2*1) 

170 WKfrAP*J) = A1 ♦ A2 

IAS*J) = A1 - A2 
IBP+J) = B1 ♦ R2 
TBN*J) a 07 — 32 
= A1 * A2 ♦ A3 

a 31 ♦ B2 ♦ B 3 

175 CONTINUF 

AJKOMjf = DCHPLX (A3 ,B3) 

K? = KO ♦ K 3 
DO 1B5 J = 1 , KH 
K1 * K1 ♦ ISP 



ORIGINAL PAGE -IS 

OF POOR QUALITY 


c 

c 


K2 
JK 
A 1 
B1 
A2 
B2 
DO 


180 

185 

190 

195 

200 

205 

210 

215 


180 
A 1 
A2 
B 1 
B2 
JK 
IP 

COHTIN 


K2 - 
J 

A0 

BO 

ZERO 
ZERO 
K 


ISP 


A (K1 ♦ 1 ) 
A(K2 + 1) 
COHTINUE 


A1 
= A2 
a B1 
= B2 
= JK 


1.KH 

♦ WK (IAP+K) 
UK (TAB+K 
BK (TBP+K 
BK (IBB+K 
J 

JK 


♦ 

f 
*♦ 
.GE. 


* 

* 

* 


BK 

BK 

BK 

BK 


ICK+JK) 

ISK+JK) 

ICK+JK) 

ISK+JK) 


KP) 


= JK - KP 


= DCBPLX (A1-B2,B1+A2) 
DCMPLX (A1+B2 f B1-A2) 


IP (KO .GT. KB) GO TO 160 
JP a KP 

IP ( I .GE. MB ) GO TO 195 
I - I ♦ 1 
GO TO 55 
I a KB 


LI 

KB 

IP 

JJ 

IP 

I 

JJ 

GO 


a 0 

= IHKfID+I-1) ♦ KB 
(KB ,<JE. KN> GO TO 215 
= IWKJIC+I-2) ♦ JJ 
(JJ .IT. IBK (IC+I-1)) 

a JJ - IBK(IC*I) 

TO 200 

(I .HE. Hfl) GO TO -210 
a Cl 

a CB * Cl - 
a SB * C2 ♦ 

TO 70 


GO TO 205 


SB 

CM 


SI 

SI 


^IWKpL+U 
» 1 


220 

225 

230 


235 

240 


Cl 
SI 
GO 
IP 
GO 
I 

JA a KT - 1 
KA = JA ♦ 1 
IF(JA.LT. 1) GO TO 
DO 220 II a 1, JA 
J a KA - II 
IWK(J+1) = IBK (J+1) - 
I a IHK (J* 1) ♦ I 

COHTIHOE 


.EQ. 1) I a I ♦ 1 


225 


IP 

J 

I 

KB 

K2 

K3 

JJ 

JK 

KO 


(KT 
i 1 
■ 0 
a 0 


.LB. 0) GO TO 27 0 


THE RESULT IS HOB 
NORHAL ORDER. 


PERBUTED TO 


IBK (ID+J) ♦ KB 
K2 

IBK (IC +J-1 ) 

JJ 

KB ♦ JJ 
ISP = IBK (IC+ J) - JJ 
K a KO ♦ JJ 
7A4 a A (K0+ 1 
A (K0+1| 


- 

A (K0< 

A (K2< 

KO = KO 
K2 = K2 
IP (KO 
KO = KO 
K 2 = K2 


(K0+ 1) 
a A(K2 + 
= ZA4 

♦ 1 
♦ 1 

, LT. K) 

♦ ISP 

♦ ISP 


GO TO 240 


156 


no nno 


157 


I? (KO .LT. K3) GO TO 235 
IP (KO .GB. K3 ♦ ISP) GO TO 2&5 
KO = KO - IffK(ID*J) ♦ JJ 
GO TO 235 

2*5 K3 = Iff K (ID*J) ♦ K3 

IF (K3 - KB .GE. IffK (ID*J- 1) ) GO TO 250 

K2 = K3 * JK 

JK = JK * JJ 

KO = K3 - IffK (ID*J) * JK 

GO TO 235 

250 IF (J .GS. KT) GO TO 260 
K a IffK (J*1) ♦ I 
J a J ♦ 1 

255 T * I ♦ 1 

- _ j 

K) GO TO 255 


0) GO TO 265 




IffK (ILL+I) 

IF JI .LT. 

GO TO 230 
260 KB = JC3 

IF (I .LE. ^ 

J = IffK (ILL *1) 

I a i - 1 

GO TO 230 
265 IF (KB .GE. K) 

J a 1 
GO TO 230 

270 JK = IffK(IC*KTL 
ISP a IffK (~~ 

H a H - 

KB = ISP/JK-2 
IF (KT .GB. B-1 ) 
ITA a ILL*K3*1 
ITB a ITA*JK 
IDN1 a ID-1 
IKT a KT*1 
IN a H*1 

DO 275 J = IKT, IN 
IHK(IDN1*J) a 
275 CONTINUE 
JJ a 0 

DO 290 J a 1,KB 
K 2 KT 

JJ = IHK (ID* R* 1) 
IF (JJ .IT. IffK 
JJ = JJ - IffK (I 
K = K ♦ 1 
GO TO 280 


GO TO 270 


GO TO 9005 


IffK (IDH1 ♦ J) /JK 


280 


£ 8 ,*”' 


GO TO 285 


285 


IF (JJ 
290 CONTINUE 


IffK (ILL* J) 
"■ EQ. 


= JJ 

J) IffK (ILL* J) = -J 


DO 300 J a 1,KB 


DETERMINE THE PER NOTATION CYCLES 
OF LENGTH GREATER THAN OR EQOAL 
TO TffO. 


295 


IF (IffK (TLl* J) .LB. 0) GO TO 300 
K2 = J 

K2 a TABS (IffK (ILL*K2)) 

IF (K2 .EQ. J) GO TO 300 
IffK (ILL* K*) = -I WK (ILL* K2) 

^0 295 


GO 


300 CONTINUE 


I a 0 

J = 0 
KB = 0 
KM = N 

3 05 J a J ♦ 1 

IF (IffK (I LT.*J) 


REORDER A FOLLOWING THE 
PERMUTATION : YCLES 


•LT. 0) GO TO 305 
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K = 'IWK (ILL*J) 

KO = JK * K ♦ KB 
3 10 ZA4 = A (K0 + I+ 1) 

WK(ITA*T) = A 4 
WK (ITB^Ij = B4 
I = I ♦ 1 

IP (T .LT. JK) GO TO 310 
I = 0 

315 K = -IWK (ILL+K) 

JJ = KO 

KO = JK * K ♦ KB 
320 A (JJ 1) = A (K 0 I ♦ 1 ) 

I » I ♦ 1 

IF (I .LT. JK) GO TO 320 
1 = 0 

IP (K . NE. J) GO TO 315 
325 A(KO*I*1) = DCMPLX (UK (ITA+I) , WK (ITB*I) ) 

IP~(I tlT. JK) GO TO 325 
1=0 

IP {J .LT. K2) GO TO 305 
J = 0 

KB = KB ♦ ISP 
IP (KB .LT. KN) GO TO 305 
9005 RETURN 
END 

IHSL 300TINE NAME - PPT2C 


COMPUTER - IBH/DOUBLE 

LATEST REVISION - JANUARY 1, 1978 

PORPOSE - COMPOTE THE PAST POOR IE R TRANSPORN CP A 

COMPLEX VALOED SEQUENCE OP LENGTH EQUAL TO 
A POWER TWO 

USAGE - CALL PFT2C (A,H,IWK) 

ARGUMENTS A r COMPLEX VECTOR OP LENGTH N, WHERE N=2**M. 

ON INPUT A CONTAINS THE COMPLEX VALUED 
SEQUENCE TO BE TRANSFORMED. 

ON OUTPUT A IS REPLACED BY THE 
FOURIER TRANSFORM. 

H - INPUT EXPONENT TO WHICH 2 IS RAISED TO 

PRODUCE THE NUMBER OF DATA POINTS, N 
(I.E. N = 2**H). 

IWK - WORK VECTOR OP LENGTH H*1. 

PRECISION/HARDW ARE - SINGLE AND DOUBLE/H32 

- SIN GL E/H36 , H4 8, H50 

REQD. IMSL ROUTINES - NONE REQUIRED 

NOTATION - INFORMATION ON SPECIAL NOTATION AND 

CONVENTIONS IS AVAILABLE IN THE MANUAL 
INTRODUCTION OR THROUGH IHSL ROUTINE UHELP 


REMARKS 1. 


2. 


PFT2C COMPUTES THE POURIER TRANSPOSE, X, ACCORDING 
TO THE FOLLOWING PORMULA; 


X(EO) = 
FOR K=0, 1 


SUM FROM J = 0 TO N-1 OP 
A (J* 1) *CEXP { (0.0. (2 . 0*PI *J*K) /N) ) 
.,N-1 AND PT-S.1'415... 


NOTE THAT X OVERWRITES A ON OUTPUT. 
PFT7C CAN BE USED TO COMPUTE 
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X(K+1) = (1/N)*SUM FROM J = 0 TO N-1 OP 

A (J*1J *CEXP(( 0. 0. (-2. 0*PI*J*K) /N) ) 

FOR K=0,1,...,N-1 AND PI=3. 1415... 

BY PERFORMING THE FOLLOWING STEPS; 

DO 10 1=1, N 

A (I) = CONJG (A (I) ) 

10 CONTINUE 

CALL PFT2C (A,M,IWK) 

DO 20 1=1, N 

A (I) = CONJG (A (I))/N 
20 CONTINUE 

- 1973 BY IHSL, INC. ALL RIGHTS RESERVED. 

- IHSL WARRANTS ONLY THAT IHSL TESTING HAS BEEN 

APPLIED TO THIS CODE. NO OTHER WARRANTY, 
EXPRESSED OR IMPLIED, IS APPLICABLE. 


SUBROUTINE FFT2C (A,H,IVK) 

- SPECIFICATIONS POR ARGUMENTS 

H, IWK(1) 

A ( 1) 

SPECIFICATIONS FOR LOCAL VARIABLES 

I, ISP,J.JJ, JSP,K,K0,K1,K2,S3,KB,KN,HK, MH,BP,N, 

& ‘if «•*»-» 1 > 

ZERO,6ne:ZO (21 .Z 1 T2) ,Z2 (2) ,Z3 (2) 

zao.4ai,4a2,zX5.aK2 


IHTEGER 
COMPLEX* 16 

( INTEGER 

DOUBLE PRECISION 


CQHPLEX*16 

EQUIVALENCE 


(ZAO ,Z0 (1) ) , (ZA i , Z 1 (1) ) , <ZA2.Z2(1) ) f 


DATA 


B3,Z3 

M *1 A 


m 


7811 865475D0/, 


DATA 


MP = M*1 
N = 2**H 
TWK(1) = 1 

55 : iVi 2) * 2 


QA70' 

SK/. 382683 4323 6 5 08 9B DO/, 

CK/. 92 3879 53 25 1 1 2868D0/. 

TWOPI/6. 28318530717958650/ 

ZERO/0. ODOAONE/1. 0D0/ 

Sg=SJPT2^2,SK=SU (PI/8) ,CK=COS (PI/8) 

FIRST EXECUTABLE STATEMENT 


DO 5 1=2, MP 

= i»MI-i)*xwk(i-1) 

CONTINUE 
RAD * TWOPI/N 
BK = M - 4 
KB = 1 

IF (MM .EQ. M) GO TO 15 
K2 = KN 

KO = IWE(Bn*1) ♦ KB 
K2 = K2 - 1 
KO = KO - 1 
AK2 = A (K2) 

A (K21 = A(K0) - AK2 
AjKOi = AjKO) *■ AK2 
... I? (KO .GT. KB) GO TO 10 
15 Cl = ONE 


INITIALIZE WORK VECTOR 


10 
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c 


c 

c 


SI = 2ERO 
JJ = 0 
It # M • 1 
.7=4 

IP (K . GE. 7) GO TO 30 
GO TO 70 

20 IP (IWK (J) .GT. JJ) GO TO 25 
JJ = JJ - IWK(J) 

IP (IUK(J) .GT. JJ) GO TO 25 
JJ = JJ - I OK (J) 

J = J - 1 
K = K «• 2 
GO TO 20 

25 JJ = IWK (J) ♦ JJ 
.7=4 

30 ISP = IVK(K) 

IP (JJ .EQ. 0) GO TO 40 



C2 

3 

J.T * 

ISP * RAD 



Cl 

= 

DC OS 

(C 2) 



SI 

3 

DSIN 

(C2) 


35 

C2 

3 

Cl * 

Cl - SI * 

SI 


S2 

3 

Cl * 

(SI ♦ SI) 



C3 

3 

C2 * 

Cl - S2 * 

SI 


S3 

3 

C2 • 

SI ♦ S2 • 

Cl 

40 

JSP 

= ISP 

♦ KB 



RESET TRIGONOMET1 IC PARAMETERS 


DETERMINE POORIER COEFFICIENTS 
IN GROUPS OP 4 


DO 50 I=1,ISP 
KO a JSP - I 


45 


K1 = KO 
K2 = K 1 
K3 = K2 
ZAO » A 
ZA1 = A 
ZA2 
Z A3 


♦ ISP 

♦ ISP 

♦ ISP 
KO) 

.. . K1 

A JK2 
A(K3 


IP (SI .EQ 
»P = A1 


TEMI 
A 1 = 
B 1 


ZERO) GO TO 45 


A1 * 
TEMP 


TEMP = A2 
A2 = A2 * . 

B2 = TEMP * S2 ♦ 


Cl - B 1 
* Si ♦ 

C2 - B2 


TEMP = A3 


A3 
B3 
TEMP 


A3 
TEMP 
AO 


* C3 - B3 


A2 = AO - 
AO = TEMP 
TEMP = A 1 
A3 = A1 - 
A 1 = TEMP 
TEMP = BO 
B2 = BO - 
BO = TEMP 
TEMP » B1 ♦ D3 
B3 = B1 - B3 
31 = TEMP 
DCMPLX 


» S3 

♦ A2 
A2 

♦ A3 
A3 

♦ B2 
B2 


* SI 

B 1 * Cl 

* S2 
B2 * C2 

* S3 
B3 * C3 


A (KO 
A K 1 
A(K2 
A K3, 

50 CONTINUE 
IP (K .LE. 
K = K - 2 
GO TO 30 


DCMPLT 

OCMPT.X 

DCMPI.X 


AO 

AO 

a: 

A2 


♦A1,B0*B1 
-A1 ,B0-B1 
- 83* B 2*A3 
♦B3,B2-A3 


1) GO 70 55 
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55 KB = K3 ♦ ISP 


I? (KH .LB. KB) GO TO 70 
I? (J . NS. t) GO TO 6 0 
K = 3 
J » HK 
GO TO 20 
J = J - 1 
C2 = Cl 

17. (J .NS. 2) GO TO 65 
Cl = Cl * CK ♦ SI * SK 
SI = SI * CK - C2 * SK 
GO TO 35 


65 Cl a (Cl - SI) 
SI = (C2 + SI) 
GO TO 35 
70 CONTINUE 


IP(H .LE. 1) GO TO 9005 
HP = (!♦ 1 
JJ = 1 

INK( 1) a 1 
DO 75 I = 2. BP 

INK (I) = INK (1-1) * 2 
75 CONTINUE 

Nft » INK (BP-2) 

IP (S .GT. 2) N8 = INK (BP-3) 
N2 = INK (HP-1) 

LB a H2 

NN a INK ( HP) ♦ 1 
HP a HP-4 

J = 2 

80 JK a jj ♦ H2 
AR2 = A (J) 

A (J) = A(.TK) 

A(JK) - AK2 
J = J*1 

IP (JJ .GT. N4) GO TO 85 
JJ = JJ ♦ N4 
GO TO 105 
85 JJ = jj - Ntt 

IP (JJ .GT. H8) GO TO 90 
JJ a JJ ♦ H8 

GO TO 105 
90 JJ » JJ - H8 

K a BP 

95 IP (INK (K) .GE. JJ) GO TO ICO 
JJ a jj - IWK(R) 

K a K - 1 
GO TO 95 

100 JJ a IWK(K) *■ JJ 
105 IP (JJ .LSI J) GO TO 110 
K = NN - J 

jk = m - jj 
AK2 = A (J) 

A (J) = A ( JJ) 

Ajjj) = AK2 
AK2 = A (K) 

A (K) = A ( JK) 

A (JK) = AK2 
110 J = J ♦ 1 


CHECK POR CQHPLETIOH OP FINAL 
ITERATION 


PERMUTE THE COMPLEX VECTOR IN 
REVERSE BINASI ORDER TO NORHAL 
ORDER 


INITIALIZE NORK VECTOR 


DETERMINE INDICES AND SNITCH A 


CYCLE REPEATED UtJTTL LIMITING NUMBER 
O? CHANGES IS ACHIEVED 
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IP (J . L E. LMJ GO TO 30 


9005 RETURN 
END 

IMSL ROUTINE NAME 


- 7FT3D 


COMPUTER 
LATEST REVISION 
PURPOSE 

OSAGE 

ARGUMENTS A 


IA1 

IA2 


N1 

N2 

N3 

I JOB 


INK 

RWK 

CWK 


PRECISION/HARDN ARE 

REQD. I MSL ROUTINES 
NOTATION 


- IBM/DOUBLE 

- JUNE 1, 1980 

- COMPUTE TEE PAST FOURIER TRANSFORM OP 

A COMPLEX VALUED 1,2 OR 3 DIMENSIONAL 
ARRAY -- 

- CALL PFT3D (A,IA1,IA2,N1,N2,N3,IJOB,IWK,RWK, 

CWK) 

- COMPLEX ARRAY. A MAY BE A THREE 

DIMENSIONAL ARRAY OP DIMENSION N1 BY N2 
BY N3, A TWO DIMENSIONAL ARRAY OP 
DIMENSION N1 BY N2, OR A VECTOR OF 
LENGTH N1. ON INPUT A CONTAINS THE 
ARRAY TO BE TRANSFORMED. ON OUTPUT 
A IS REPLACED BY THE FOURIER OR 
INVERSE FOURIER TRANSPOSE (DEPENDING OH 
THE VALUE OF INPUT PARAMETER IJOB) . 

- FIRST DIMENSION OF THE ARRAY A EXACTLY 

AS SPECIFIED IN THE DIMENSION STATEMENT 
IN THE CALLING PROGRAM. (INPUT) 

- SECOND DIMENSION OF THE ARRAY A EXACTLY 

AS SPECIFIED IN THE DIMENSION STATEMENT 
IN THE CALLING PROGRAM. (INPUT) 

- LIMITS ON THE FIRST, SECOND, AND THIRD 

SUBSCRIPTS OF THE ARRAY A, RESPECTIVELY. 
(INPUT) 

- INPUT OPTION PARAMETER. 

IF IJOB IS POSITIVE, THE PAST FOURIER 
TRANSFORM OP A IS TO BE CALCULATED. 

IF IJOB IS NEGATIVE, THE INVERSE 

FAST FOURIER TRANSFORM OF A IS TO BE 
CALCULATED. 

- INTEGER WORK VECTOR OF LENGTH 

6*MAX(N1,N2,N3) *150. 

- REAL WORK VECTOR OF LENGT1 

6*MAX(N1.N2,N3) +150. 

- COMPLEX WOftK VECTOR OF LENGTH 

MAX (N2,N3) . 

- SINGLE AND DOUBLE/H32 

- SINGLE/H32, H48, H60 

- PFTCC 

- INFORMATION ON SPECIAL NOTATION AND 

CONVENTIONS IS AVAILABLE IN THE MANUAL 
INTRODUCTION OR THROUGH IHSL ROUTINE UHELP 


REMARKS 1. IF IJOB IS POSITIVE, FFT3D' CALCULATES THE FOURIER 
TRANSFORM, X, ACCORDING 70 THE FOLLOWING FORMULA 

X(I*1, J* 1, K* 1) =TRIPLE SOM OF A (L* 1 , N ♦ 1 , N* 1) * 

Ex Pf 2*PI*SQRT (-1) * (I*L/N1+J*M/N2 *K*N/n3) ) 

WITH L=0. . . N1-1 , N=0...N2-1, N=0. - . N3-1 
AND PI=3. 1415... 


IF IJOB IS NEGATIVE. FFT3D CALCULATES THE INVERSE 
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POURIEP TRANSFORM, X, ACCORDING TO THE FOLLOWING 
FORMULA 

X (1+1, J*1, K+ 1) =1/(N1*N2*N3) *TRIPLE SOW OF 

EXP(-2*PI*SQRT(-1) *(I*L/N1+J*H/N2+K*N/N3)) 

WITH L=0...81-t, h=0. . .N2-1 , N=0...N3-1 
AND PI=3. 1415. .. 

NOTE THAT X OVERWRITES A ON OOTPOT. 

IF A IS A TWO DIMENSIONAL ARRAY, SET N3 = 1. 

IF A TS A ONE DIMENSIONAL ARRAY (VECTOR) , 

SET IA2 = N2 * N3 = 1. 

- 1980 BY IHSL, INC. ALL RIGHTS RESERVED. 

- IHSL WARRANTS ONLY THAT IHSL TESTING HAS BEEN 

APPLIED TO THIS CODE. NO OTHER WARRANTY, 
EXPRESSED OR IMPLIED, IS APPLICABLE. 


C 

c 


SUBROUTINE FFT3D 
INTEGER 

DOUBLE PRECISION 
COMPLEX* 16 

INTEGER 

DOUBLE PRECISION 
COMPLEX*16 


(A, IA 1, IA2, SI ,N2 , H3 . I JOB, IWK. RWK ,CWR) 

SPECIFICATIONS FOR ARGUMENTS 
IA1,IA2,N1,N2,N3,IJOB,IWK (1) 


RWK(I) 

(ill. 


A 


C123 

IP (IJOB.GT.O) GO TO 10 

DO 5 1=1, N1 
DO 5 J= ij N2 
DO 5 K= 1 * N3 


I,J § K,L,H,N 


IA2,N3^CWKJ1^ 


SPECIFICATIONS FOR LOCAL VARIABLES 


FIRST EXECUTABLE. STATEHENT 
INVERSE TRANSFORM 


5 CONTiN6E ,Kl 


DCONJG (A(I, J,K) ) 

TRANSFORM THIRD SUBSCRIPT 


10 DO 25 L=1, HI 
DO 25 H=1,N2 
DO 15 N=1,N3 

CWR(N) = A (L, H, N) 

15 CONTINUE 

CALL FPTCC (CWK , N3 , IWK, RWK) 

DO 20 F=1,R3 

A(L,H,K) = CWK (K) 

20 CONTINUE 

25 CONTINUE 

TRANSFORM SECOND SUBSCRIPT 

DO AO L=1 , N 1 
DO 40 K=1,N3 
DO 30 H=1,N2 

_ A CWKjf.H) = A (L, H, K) 

30 CONTINUE 

CALL FFTCC ( CWK, N2, IWK, RWK) 

DO 35 J=1, N2 

a]L.J,K) = CWK (J) 

35 CONTINUE 
40 CONTINUE 

TRANSFORM PIRST SUBSCRIPT 

DO 45 J=1,N2 
DO 45 K=1,N3 

CALL FFTCC ( A ( 1 , J, K) ,N 1 , TW K, RWK) 

45 CONTINUE 

TP (IJOB.GT.O) GO TO 55 
R 12.1 = N1*N2*N3 



C 123 = DC3PLX (R123,0.0D0) 

DO 50 I* 1, Ml 
DO 50 J=1,N2 
DO 50 K=1,H3 

50 = »“»«U(I.J.Kn/C123 

55 RETURN 
END 



